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1.   Introduction 

This  review  is  intended  to  summarize  the  developments  in  numerical  weather 
prediction  which  have  occurred  during  approximately  the  last  three  years.   This 
roughly  covers  the  period  since  the  completion  of  two  major  texts  by  Haltiner 
(1971)  and  Phillips  (1973).   Earlier  research  will  be  referenced  only  when 
required  for  the  development  of  a  particular  topic.   An  indication  of  the  in- 
tense activity  in  this  field  can  be  seen  in  the  papers  presented  at  the  Second 
American  Meteorological  Society  Conference  on  Numerical  Prediction  in  Monterey, 
California  during  1-4  October  1973,  for  which  the  abstracts  have  been  published 
in  the  July,  1973  Bulletin  of  the  American  Meteorological  Society. 

Since  the  first  application  of  numerical  integration  techniques  to  large 
scale  weather  prediction  [Charney,  Fjortoft,  and  von  Neumann  (1950)],  these 
techniques  have  been  applied  to  a  large  variety  of  atmospheric  and  oceanic  pre- 
diction problems.  These  include  studies  of  the  general  circulation  of  the 
atmosphere  and  the  oceans,  tropical  storms,  fronts,  sea  breezes,  cumulus  con- 
vection, atmospheric  pollution,  the  planetary  boundary  layer  and  turbulence. 
This  review  will  focus  on  the  developments  in  large  scale  weather  prediction, 
although  many  of  the  results  have  applications  to  other  problems.  The  para- 
meterization of  smaller  scale  phenomena  will  not  be  considered  in  this  review 
because  it  deserves  a  separate  comprehensive  treatment. 

Numerical  integration  of  the  atmospheric  equations  as  an  initial  value 
problem  is  the  primary  basis  for  the  prediction  of  synoptic-scale  disturbances 
for  periods  between  12  hours  and  perhaps  five  days  and,  in  addition,  to  some 
extent  for  smaller  scales  and  much  longer  periods.  The  sources  of  error  in 
such  predictions  are  a  consequence  of  (a)  gaps  and  errors  in  the  data  which 


make  up  the  initial  state,  (b)  limitations  in  the  objective  analysis-initiali- 
zation schemes  which  are  applied  to  the  data,  (c)  truncation  errors  in  numeri- 
cal integration  schemes,  (d)  incomplete  representation  of  the  many  complicated 
dynamical  processes  at  work  in  the  atmosphere  and  finally,  limitations  imposed 
by  the  predictability  of  the  atmosphere. 

2.   Objective  Analysis  and  Initialization 

A  standard  method  of  objective  analysis  which  has  been  in  use  for  nearly 
15  years  or  so  begins  with  an  initial  guess  based  on  a  prognosis  of  the  mass 
field  from  the  previous  analysis  or,  perhaps  simply  a  persistence  forecast. 
The  initial  gridpoint  values  are  then  modified  with  actual  observations 
weighted  according  to  the  distance  between  the  observation  and  the  gridpoint. 
Winds  are  obtained  geostrophically  or  by  means  of  a  balance  equation.  With 
some  care  to  insure  vertical  consistency  and  to  avoid  superadiabatic  lapse 
rates,  this  type  of  objective  analysis  is  quite  adequate  for  initializing  the 
vorticity  (filtered)  type  of  prediction  models.  However,  the  primitive  equa- 
tion or  P.E.  models  introduced  operationally  about  a  half  dozen  years  ago, 
are  much  more  sensitive  to  the  initial  conditions  than  are  the  filtered  models 
from  which  gravity  waves  have  been  eliminated.  As  a  consequence,  the  lack 
of  proper  balance  between  the  pressure  and  wind  fields  can  give  rise  to  quite 
large,  spurious  inert ial-gravity  waves  with  periods  ranging  from  a  few  time 
increments   At  to  a  day  or  more.  Through  the  natural  geos trophic  adjustment 
process,  which  is  inherent  in  the  P.E.  models,  these  spurious  gravity  modes 
are  eventually  dispersed  and  suppressed  to  acceptable  magnitudes,  perhaps 
comparable  to  those  observed  in  nature.  To  avoid  the  large  pressure  fluctua- 
tions associated  with  these  spurious  gravity  waves,  greater  care  should  be 


taken  to  bring  the  objective  analysis  to  a  state  where  the  wind  and  pressure 
fields  are  more  nearly  in  proper  balance  before  beginning  the  forecast,  a 
procedure  which  is  usually  referred  to  as  initialization. 

The  first  efforts  in  this  direction  simply  followed  the  practice  in  fil- 
tered models  by  using  winds  from  a  stream  function  obtained  by  solution  of 
the  balance  equation  with  the  objectively-analyzed  geopotential  fields  on 
constant  pressure  surfaces.   Since  the  balance  equation  allows  for  curved 
motion  to  some  extent,  the  resulting  winds  are  in  better  balance  with  the 
pressure  force  than  geostrophic  winds.   Nevertheless,  as  shown  by  Phillips 
(1960),  an  appropriate  divergent  wind  component  is  necessary  if  the  inertial- 
gravity  waves  are  to  be  excluded  or  negligible.   Consequently,  the  rotational 
winds  obtained  by  solution  of  the  conventional  balance  equation  are  not 
sufficient  to  prevent  the  unwanted  gravity  modes  in  a  primitive  equation 
model,  although  they  are  clearly  better  than  geostrophic  winds  or  unmodified 
observed  winds. 

Recently,  Sundqvist  (1973),  on  the  basis  of  limited  tests,  suggested  that 
the  solution  of  a  balance  equation  with  zero  mass  divergence  (V   •  TT\V  =  0) 
on  sigma  surfaces  would  lead  to  initial  winds  which  give  less  gravitational 
noise  than  the  winds  from  the  conventional  balance  equation  on  pressure  sur- 
faces followed  by  interpolation  to  sigma  surfaces. 

In  another  approach,  Temperaton  (1973)  verified  that  gravity  oscillations 
were  indeed  present  after  initialization  with  the  exact  rotational  part  of 
the  wind.   The  latter  was  obtained  from  an  idealized  experiment  where  "con- 
trol" data  was  first  generated  by  first  running  a  numerical  prediction  model 
for  two  days  using  the  Euler-backward  scheme  to  dampen  the  high  frequency 


modes  leaving  essentially  only  the  slowly  evolving  meteorological  modes. 
Then  the  geopotential  field  was  treated  as  the  observed  field  and  a  series 
of  experiments  were  run  for  an  additional  day  with  different  initializa- 
tion methods.   In  particular  the  rotational  wind  component  was  extracted 
from  the  control  field  by  solving  the  Poisson  equation,  v  f  =  k  •  Vx  V  , 
for  the  stream  function  \J/  .   Initializing  with  this  wind  resulted  in 
gravity  noise  comparable  to  that  resulting  from  winds  obtained  from  solu- 
tion of  the  balance  equation. 

The  divergent  component  of  the  wind  implied  by  the  quasi-geos trophic 
theory  can  be  obtained  by  solving  the  corresponding  co-equation  and  then 
using  the  continuity  equation  to  obtain  the  divergent  wind  potential  X 
as  follows: 


■I 


X=-^T  (1) 


This  is  a  rather  lengthy  task  however,  and  is  probably  less  effective  and 
maybe  as  time  consuming  as  simply  using  the  model  equations  to  achieve 
the  desired  balance.   In  fact,  Houghton,  Baumhefner  and  Washington  (1971) 
made  a  study  of  initialization  for  the  NCAR  global  model  and  found  that 
inclusion  of  the  vertical  velocity  from  an  co-equation  and  the  corresponding 
divergent  wind  component  in  the  initial  winds  did  not  significantly  reduce 
the  trauma  of  inertial  gravity  oscillations. 

At  the  National  Meteorological  Center  (NMC)  solution  of  the  balance 
equation  was  abandoned  several  years  ago  in  favor  of  extracting  the  rota- 
tional wind  component  directly  from  objectively  analyzed  winds.   In  addi- 
tion, the  12-hour  forecast  divergent  wind  component  is  utilized.  The 


resulting  total  initial  wind  field  gives  better  agreement  with  the  observa- 
tions without  significantly  increasing  the  noise  level  in  the  prediction 
model. 

Although  they  are  usually  treated  as  separate  steps,  objective  analysis 
and  initialization  have  the  same  goal,  namely,  to  provide  an  accurate, 
properly  balanced  state  from  which  the  hydrodynamical  equations  can  be 
numerically  integrated  forward  in  time.   Obviously,  the  two  steps  are  not 
dynamically  independent  and  need  not  be  separated,  which  is  the  view  taken 
by  a  number  of  investigators.  A  specific  example  is  Sasaki's  (1958)  ini- 
tialization by  variational  methods  which  incorporates  dynamical  constraints 
beyond  the  usual  geostrophic  wind  law  in  the  objective  analysis.  These 
constraints  have  included  the  use  of  a  generalized  wind  equation,  suppres- 
sion of  high-frequency  oscillations,  as  well  as  minimizing  RMS  differences 
between  observations  and  a  first  guess  field. 

Lewis  (1972)  applied  Sasaki's  technique  to  develop  an  operational  analy- 
sis of  wind  and  temperature  for  the  global  band  40  S  to  60  N  from  the 
surface  to  250  mb.  Disregarding  vertical  advection,  Lewis  derived  the 
following  generalized  thermal  wind  equation 

St  53  "  "  55  (v  '  W)  -  m   x  55  "  RVI  s  F         (2> 

The  variational  formulation  calls  for  the  minimization  of  the  integral 
I  =  fff  [a(v  -  V)2  +  P(T  -  T)2  +  a(Fx2  +  Fy2)]       (3) 

where  a,  P  are  precision  moduli  for  the  wind  and  temperature  based  on 
variances  and  a  is  a  dynamical  weight  factor.  Solution  of  the  Euler 


equations  associated  with  the  above  integral  provide  the  new  wind  and  tem- 
perature fields.  At  the  Monterey  NWP  meeting  in  October  1973,  Sasaki 
applied  a  similar  procedure  to  the  problem  of  initialization  by  utilizing 
a  dynamic  constraint  on  the  kinetic  energy  which  reduced  the  gravity  wave 
noise  considerably. 

The  variational  principle  has  also  been  used  by  Flattery  (1967)  at  NMC 
in  making  objective  analyses  by  fitting  Hough  functions  (the  eigenfunctions 
of  Laplace's  tidal  equation)  to  the  observed  data  by  a  least  squares  method. 
The  results  have  been  used  in  a  global  spectral  model. 

Another  common  method  of  objective  analysis  interpolates  between  ob- 
servations to  gridpoints  by  requiring  the  gridpoint  values  to  satisfy  a 
Poisson  equation.   Starting  with  an  initial  guess,  the  gridpoint  values  are 
then  modified  by  the  use  of  observed  values  which  are  treated  as  internal 
boundary  points.  This  scheme  was  evaluated  by  Leary  and  Thompson  (1973) 
for  accuracy  with  known  spherical  harmonics.  Wavenumber  2  was  reasonably 
well  represented  with  87^  of  its  amplitude  squared  appearing  in  the  analysis 
field,  while  only  13$  of  the  input  amplitude  squared  of  wavenumber  12  sur- 
vived the  analysis.  Aside  from  demonstrating  a  rather  severe  deficiency 
of  this  analysis  scheme,  the  study  suggested  that  there  is  a  less  steep  drop- 
off in  the  kinetic  energy  of  the  wind  field,  perhaps  a  "-2"  rather  than  the 
"-3"  power  law  obtained  from  spectra  of  such  objectively  analyzed  data. 

Several  other  approaches  have  been  utilized  to  suppress  the  inertial- 
gravity  oscillations  usually  present  in  the  early  stages  of  a  P.E.  forecast. 
One  procedure  consists  of  integrating  forward  and  backward  about  the  initial 
time  starting  with  the  objectively  analyzed  data  and  periodically  restoring 


either  the  mass  or  the  wind  field  at  the  initial  time,  usually  the  mass  field 
in  middle  latitudes  and  perhaps  the  wind  field  in  tropical  latitudes.   For 
this  purpose  it  would  appear  advantageous  to  use  an  integration  method  which 
selectively  damps  high-frequency  oscillations  such  as  the  Euler-backward 
scheme  introduced  by  Matsuno  (1966).   As  applied  to  a  simple  advective  equation 

The  Euler  backward  scheme  will  damp  short  waves  about  10$  with  each  time  step, 
while  long  meteorological  type  waves  will  be  negligibly  damped.   However, 
internal  gravity  waves  of  low  frequency  may  also  dampen  rather  slowly.  Thus 
considerable  computing  time  may  be  required  by  this  two-step  scheme  to  accom- 
plish the  desired  goal  of  suppressing  the  spurious  gravitational  noise,  perhaps 
the  equivalent  of  a  24-  or  48-hour  forecast  which  is  operationally  undesirable. 

The  presence  of  the  gravity  modes  is  usually  reflected  in  the  horizontal 
divergence,  the  vertical  velocity  and  the  surface  pressure  tendency.  Recog- 
nizing that  the  divergent  part  of  the  wind  is  intimately  related  to  the 
propagation  of  the  gravity  oscillations,  Talagrand  (1972)  suggested  a  special 
viscosity  term  to  dampen  the  divergent  wind  component  as  follows: 

du     du  dD     br        A 

^  +  U^+ -^+V^=° 

(5) 

Here  D  is  the  horizontal  velocity  divergence  and  Q    ,  the  relative  vorticity. 
Now  when  the  divergence  and  vorticity  equations  are  formed,  one  obtains 

t  + -  H  ^D  =  0 

(6) 
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It  is  clear  that  with  V  =  0  and  (_i  ^  0  ,  divergence  but  not  vorticity  will 
be  suppressed.   McPherson  (1973)  of  NMC  evaluated  this  novel  concept  with  an 
8-layer  global  primitive  equation  model  using  Shuman's  semi-momentum  differ- 
encing. With  (j  =  0  ,  oscillations  primarily  in  the  form  of  external  gravity 
modes  appeared  in  middle  latitudes  as  follows: 

AMPLITUDE  PERIOD 

5-10  meters  3-4  hours 

20  -  25  meters  6-7  hours 

20  -  60  meters  10  hours 

These  waves  were  largely  eliminated  with  a  viscosity  coefficient  of  |j  =  2.5 

8  2 
x  10  m  /sec  in  about  7  hours  time.   Several  unfortunate  side  effects 

occurred; however,  there  was  abnormally  high  surface  pressure  at  the  end  of 

12  hours  in  the  vicinity  of  mountains  together  with  a  downstream  trough. 

7  2 
Also  values  of  p.  in  excess  of  10  m  /sec  eliminated  all  precipitation. 

Clearly  further  evaluation  is  necessary.  McPherson  speculated  that  the 
slope  of  the  sigma  surfaces  in  the  vicinity  of  mountains  may  be  responsible. 
We  would  not  find  this  surprising  since  two-dimensional  diffusion  on  steep- 
sloping  sigma  surfaces  near  mountains  may  involve  large  vertical  wind  shear 
and  temperature  differencing,  with  undesirable  consequences.   Perhaps  the 
diffusion  technique  would  be  more  successful  if  carried  out  on  quasi- 
horizontal  surfaces,  or  in  such  a  way  as  to  avoid  changing  the  selective 
vertical  diffusion. 

At  this  point  we  would  like  to  return  to  Temperaton's  experiment  on  dy- 
namical initialization.  He  used  the  so-called  "shallow  water"  equations  in 
flux  form  with  spherical  coordinates.   In  previous  experiments  several  years 
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ago  Nitta  and  Hovermale  (1967)  had  initialized  by  integrating  forward  N  time 
steps,  backward  2N  steps,  and  then  returning  to  the  initial  time  with  N  steps 
forward;  however,  N  was  taken  to  be  just  1.   The  Euler-backward  scheme  gave 
somewhat  faster  convergence  than  the  leapfrog  scheme.  Mesinger  (1972)  recom- 
mended a  larger  amplitude  NAt  and  partial  restoration  of  mass  more  frequently. 
Winninghoff  (1971)  integrated  a  barotropic  model  18  hours  forward  in  time  and 
back  to  the  initial  time  with  the  2-step  Euler-backward  scheme  and  obtained 
convergence  to  quite  good  initializing  winds,  but  obviously  at  great  expense 
in  computer  time  (equivalent  to  a  72-hr  forecast  with  the  leapfrog  scheme). 
Neither  friction  nor  heating  terms  were  present  during  initialization.   Fric- 
tion and  simulated  heating  after  initialization  gave  somewhat  poorer  results 
than  in  the  frictionless,  adiabatic  case.   However,  including  these  effects 
in  the  initialization  procedure  could  lead  to  undesirable  results  depending  on 
the  mathematical  formulation  of  these  processes,  which  may  not  be  reversible. 

Temperaton  tried  restoring  mass  only  partially  at  the  central  time  but 
was  not  successful.   Using  a  time  amplitude  of  4At  and  2/3  restoration  of  the 
mass  field  at  every  other  time  step  gave  quite  rapid  convergence  with  respect 
to  the  total  RMS  wind  error.  However,  gravity  wave  activity,  as  measured  by 
the  RMS  divergence,  was  not  diminished.   Analysis  of  the  winds  showed  that 
the  rotational  wind  component  converged  rapidly  toward  the  correct  value  but 
the  divergent  wind  component  was  more  in  error  than  with  the  case  N  =  1.   He 
then  tried  another  approach  of  integrating  forward  N  time  steps  from  the 
initial  time,  then  N  steps  backwards  from  the  initial  time  and  finally  averag- 
ing the  end  results  as  follows: 

u(tQ)  =  |[u(tQ  +  NAt)  +  u(tQ  -  NAt)]  (7) 
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The  mass  is  now  restored  and  the  procedure  repeated  as  many  times  as  is  needed 
to  bring  about  convergence.   The  leapfrog  method  was  used  except  for  a  first 
forward  time  step.   The  averaging  tends  to  remove  high-frequency  modes  but 
leaves  low- frequency  waves  relatively  unaffected.   For  example,  the  2NAt  fre- 
quency would  be  removed  by  such  averaging  while  the  4NAt  wave  would  be  reduced 
by  about  one-half.   Iteration  cycles  of  3 At,  6 At  and  9 At  were  tried  with  the 
best  results  obtained  from  the  6At  (1  hour)  case.   The  Euler-backward  method 
was  tried  for  starting  in  each  direction  but  showed  no  improvement  over  a 
simple  forward  step.   The  principal  result  was  that  considerably  faster  conver- 
gence was  achieved  using  the  averaging  technique  than  with  the  Euler  backward 
scheme  N  =  1.  The  high  frequency  oscillations  in  the  RMS  divergence  were 
markedly  reduced  to  about  one-tenth  the  amplitude  occurring  with  the  use  of 
the  exact  rotational  wind  component  without  the  averaging  scheme.  Nevertheless, 
there  remained  an  easily  recognizable  oscillation  with  a  period  of  about  four 
hours.  Elimination  of  this  oscillation  and  further  reduction  of  the  noise 
was  achieved  by  first  adjusting  the  wind  field  while  restoring  the  mass  field 
at  the  initializing  time  and  then  repeating  the  forward  and  backward  averaging 
procedure  except  that  the  mass  field  was  averaged  while  the  newly  computed 
winds  were  restored  after  each  cycle;  thus  both  mass  and  winds  were  adjusted. 
In  all,  20  cycles,  equivalent  to  a  40-hour  forecast,  were  completed,  which  is 
again  too  long  for  operational  forecasting. 

Suppression  of  considerable  gravity  noise  has  been  achieved  at  NMC  by  us- 
ing a  running  time  average  of  the  dependent  variables  in  connection  with  the 
leapfrog  scheme  due  to  A.  Robert  (1966).   Specifically,  when  a  new  value  of 
a  predicted  variable,  say  A(t  +  At),  has  been  calculated  using  the  leapfrog 
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time  differencing,    a  weighted  average  of  the   last   three  time  values    is   cal- 
culated as    follows: 

A*(t)   =  A(t)   +  0.5v[A*(t-At)   -   2A(t)   +  A(t+At)]  (8) 

The  new  value  A*(t)  is  now  used  in  the  next  leapfrog  step  to  A(t  +  2At), 
that  is 

dA 
A(t  +  2At)  =  A*(t)  +  2At  (§f)t+At  (9) 

The  weight  factor  v  controls  the  amount  of  smoothing;   V  =  0.5  gives  the 
familiar  1-2-1  averaging.  The  simpler  averaging  function  with  A(t  -  At) 
instead  of  A*(t-At)  in  Eq.  (8)  would  completely  remove  waves  of  period  2At 
and  dampen  a  4 At  wave  by  about  one  half,  while  much  longer  periods  are  rela- 
tively unaffected.  Recall  here  that  the  spurious  computational  mode 
generated  in  the  leapfrog  scheme  has  very  nearly  a  period  of  2At  if  a  corres- 
ponding physical  mode  is  of  a  relatively  lower  frequency.   As  a  consequence, 
NMC's  averaging  procedure  is  effective  in  removing  the  computational  mode. 
Moreover,  repetition  of  the  averaging  at  every  time  step  tends  to  damp  some- 
what lower  frequencies,  including  a  considerable  part  of  the  spurious  gravity 
noise  generated  through  the  early  hours  of  a  P.E.  forecast.  Asselin  (1972) 
has  computed  the  damping  and  phase  shift  characteristics  as  a  function  of  fre- 
quency for  a  time  filter  of  this  form.   Low  frequency  meteorological  waves  show 
negligible  amplitude  change  and  phase  shift,  but  there  is  very  effective 
damping  of  computational  modes. 

In  some  recent  initialization  experiments  with  a  five-level,  global  P.E. 
model,  Haltiner  and  McCullough  (1974)  did  not  find  any  significant  reduction 
of  gravity  noise  with  initial  rotational  winds  from  a  balance  equation  solved 
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on  a-surfaces  as  compared  to  similar  winds  from  a  p-surface  balance  equation. 

This  may  be  somewhat  surprising  since  from  theoretical  considerations,  it 

might  be  expected  that  at  least  external  gravity  waves  would  be  suppressed. 

Sundqvist  assumed  zero  mass  divergence,  V  •  TTV=  0  ,  where  TT  is  surface 

pressure,  to  obtain  the  balance  equation  on  sigma  surfaces.   It  follows  from 

the  integrated  continuity  equation,  namely,  ^TT/dt  =  -\  V  •  (TTVV)da,  that 

Jo 
if  V  •  TTV  =  0  ,  the  surface  pressure  tendency  will  vanish  initially  which 

should  tend  to  suppress  external  gravity  waves.   In  any  event,  the  Haltiner- 

McCullough  initialization  experiments  with  the  (J-balanced  winds  did  not  show 

significant  improvement  in  damping  spurious  inert ial-gravity  waves. 

The  Robert  time  filter  (1966)  is  not  only  effective  in  eliminating  the 
computational  mode  associated  with  the  leapfrog  scheme,  but  also  gradually 
tends  to  dampen  frequencies  characteristic  of  the  gravitational  noise  at  the 
beginning  of  a  P.E.  forecast.  Haltiner  and  McCullough  combined  the  time 
filter  with  the  averaging  scheme  of  Temperaton  to  substantially  reduce  the 
noise  in  the  equivalent  computer  time  of  a  12-hr  forecast.  The  latter  was 
accomplished  by  twice  averaging  winds  from  a  3-hr  forecast  and  a  3-hr  hind- 
cast  and  restoring  the  mass  field  to  the  initial  values  for  each  time.  The 
winds  thus  obtained  together  with  the  original  mass  field  constituted  the 
initial  conditions  for  prediction,  which  resulted  in  a  substantial  decrease 
in  high  frequency  surface  pressure  oscillations  compared  to  the  immediate 
use  of  the  balanced  winds  for  initialization. 

The  addition  of  the  divergent  wind  component  predicted  by  the  model  from 
a  previous  analysis  to  the  present  analysis  time  should  further  aid  in  damp- 
ing the  spurious  inertial-gravity  noise  generated  by  the  introduction  of  new 
data.   This  step  costs  very  little  in  computer  time. 
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3.   Data  Assimilation 

The  advent  of  the  satellite  infrared  spectrometers  (e.g.,  SIRS),  which 
can  sound  the  atmosphere  from  space,  has  provided  a  new  and  important 
source  of  data  which  is  especially  valuable  over  sparse-data  areas.   It 
did,  however,  further  emphasize  a  problem  of  how  to  utilize  most  effectively 
upper  air  data  which  do  not  occur  at  or  very  near  the  regular  synoptic  obser- 
vation times,  a  problem  already  existing  in  connection  with  aircraft  observa- 
tions.  The  procedure  for  incorporating  off-time  reports  into  an  analysis- 
prediction  scheme  is  referred  to  as  four-dimensional  assimilation.   Clearly 
the  trauma  of  spurious  inertial  gravity  waves  associated  with  initialization 
must  be  eliminated  if  new  observational  data  is  to  be  more  or  less  continuously 
incorporated  into  a  prediction  procedure. 

Some  of  the  early  idealized  experiments  dealing  with  the  assimilation 
of  SIRS  data  by  Charney,  et  al  (1969),  Jastro  and  Halem  (1970)  and  Williamson 
and  Kasahara  (1971)  showed  that  the  simple  insertion  of  temperature  data  can 
eventually  determine  the  wind  field  and  vice  versa.   The  process  of  inserting 
some  variables  into  a  numerical  forecast  while  others  are  unchanged  has  been 
referred  to  as  updating,  which  seems  to  be  a  rather  restrictive  definition  in 
terms  of  operational  forecasting.   In  any  event,  the  degree  to  which  one  field, 
for  example,  winds,  can  be  obtained  by  updating  another,  perhaps  temperature, 
depends  on  the  predictability  of  the  model  which,  in  turn,  depends  on  natural 
instabilities,  nonlinear  exchanges,  frictional  dissipation,  etc..  Using 
model  data  to  stimulate  the  real  world,  then  introducing  errors  and  updating 
with,  quote,  "observed",  data  provides  an  indication  of  the  best  possible 
results  that  could  be  obtained  by  updating  an  operational  prediction  model 
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with  real  data.   In  actual  operational  forecasts  there  are  other  sources  of 
errors  including  the  representation  of  the  physical  processes  and  numerical 
truncation.   Naturally,  the  better  the  model,  the  less  the  predictability 
error,  but  simple  insertion  of  data  is  assuredly  not  the  best  way  to  update. 
It  would  be  quite  logical  to  modify  other  variables  by  static  or  dynamic 
balancing  and  refer  to  the  whole  procedure  as  updating  or  four-dimensional 
data  assimilation. 

Mesinger  (1972)  suggested  the  following  three  names  for  methods  of  obtain- 
ing a  proper  balance  between  the  mass  and  motion  fields: 

(1)  Static  Balancing  -  use  of  a  wind  law  such  as  the  geostrophic 
relation  or  the  balance  equation. 

(2)  Dynamic  Balancing  -  going  backward  and  forward  about  a  central  time, 
restoring  mass  or  wind  at  t  =  0  or  letting  both  vary.  Sasaki's  variational 
method  can  also  be  considered  as  a  combination  of  static  and  dynamic  balancing, 

(3)  Four-Dimensional  Balancing  in  space  and  time  -  introduction  of  data 
three-dimensionally  periodically  at  discrete  time  intervals,  including 
regular  synoptic  times. 

Hayden  (1972)  used  a  2- layer  P.E.  model  with  Shuman's  semi-momentum  dif- 
ferencing scheme  to  run  a  series  of  experiments  with  the  insertion  of  tempera- 
ture calculated  from  geopotentials  which  were  in  turn  derived  from  the  SIRS 
temperatures.   Data  were  inserted  six  times  during  a  12-hr  forecast.  Tempera- 
ture tendencies  over  a  period  of  one  orbit  were  generally  less  than  the 
observational  error;  consequently,  as  found  by  Talagrand  (1971),  too  frequent 
updating  may  be  deleterious.  Hayden  found  that  a  static  balancing,  consist- 
ing of  a  geostrophic  wind  correction,  Av,  computed  from  the  changes  in 
geopotential  A^  inferred  from  the  temperature  insertions  aided  the  geostrophic 
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adjustment  process  and  reduced  the  shock  of  updating.  He  also  attempted 
dynamic  balancing  with  N  =  2  ,  i.e.,  integrating  two  steps  forward,  4At 
backward  and  2At  forward  to  return  to  the  initial  time. 

Three  criteria  were  used  to  measure  the  successful  assimilation  of  data: 
(a)  the  updating  must  not  unduly  shock  the  model  which  he  inferred  from  a 
measure  of  the  divergent  wind  component;  (b)  does  the  model  remember  the  data 
inserted  (this  was  tested  by  reversing  the  forecast  after  12  hours  and  noting 
whether  the  hindcast  temperatures  were  nearer  to  the  inserted  values  than 
were  the  original  temperatures  that  existed  prior  to  the  updating.);  (c)  lastly 
does  the  updating  result  in  better  mass  and  wind  distributions  when  compared 
to  the  NWS  analyzed  fields?  For  this  purpose  the  12-hour  forecast  was 
followed  by  a  12-hour  hindcast.   Then  the  difference  between  the  cycled  data 
and  the  initial  data  were  correlated  with  the  difference  between  the  NWS 
analysis  and  the  initial  state;  a  positive  correlation  indicates  success. 

The  following  conclusions  were  reached: 

a.  Even  with  poor  initial  conditions,  temperature  can  be  assimilated 
without  shock  if  some  balancing  is  performed  to  aid  the  geostrophic  adjustment 
process.   Dynamic  balancing  is  sufficient  but  static  balancing  in  the  form  of 
a  simple  geostrophic  wind  correction  speeds  up  the  assimilation  process  con- 
siderably, at  least  outside  the  tropics. 

b.  Hayden  anticipates  that  under  operational  conditions,  where  the  model 
state  is  maintained  similar  to  the  observed  state,  four -dimensional  assimila- 
tion can  be  accomplished  without  time-consuming  dynamic  initialization. 

c.  Four-dimensional  data  assimilation  is  evidently  more  effective  than 
regular  objective  analysis  of  off  time  reports  that  have  been  updated  to 
regular  synoptic  times  by  Lagrangian  advection.  The  exact  details  of  the 
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surface  pressure  appear  to  be  relatively  unimportant  to  the  effectiveness 
of  the  four-dimensional  assimilation.   However,  the  SIRS-B  data  by  them- 
selves are  not  capable  of  defining  the  circulation  because  the  data  density 
is  not  sufficient  for  objective  analysis,  nor  with  the  model,  of  producing 
even  an  approximate  surface  pressure  field. 

Bengtsson  and  Gustavsson  (1971)  had  previously  found  also  that  analysis 
prior  to  insertion  of  data,  say  from  a  satellite,  leads  to  a  more  rapid 
reduction  of  error.   Talagrand  and  Miyakoda  (1971)  showed  that  a  synthesis 
technique  of  averaging  forecasts  made  from  objective  analyses  made  at 
different  times  can  reduce  the  random  errors  of  measurement  and  analysis, 
sort  of  a  Monte  Carlo  approach.  They  did  some  studies  of  inserting  data 
into  a  running  forecast  when  and  where  the  data  were  available.   If  the 
difference  between  the  predicted  values  and  the  observations  is  less  than 
or  equal  to  the  observational  error,  don't  insert  it;  there's  no  point  to 
shocking  the  system  and  uselessly  creating  spurious  inert ial-gravity  waves. 

In  summary,  experimentation  thus  far  suggests  that  one  or  several  of 
the  following  steps  are  useful  during  initialization  and  four-dimensional 
data  assimilation. 

Firstly,  at  the  regular  0000Z  and  1200Z  synoptic  times: 

a.   Perform  a  preliminary  objective  analysis  and  obtain  a  stream  func- 
tion by  solution  of  the  balance  equation,  perhaps  on  the  surfaces  of  the 
vertical  coordinate,  or  as  an  alternative  obtain  the  rotational  wind  directly 
from  objectively  analyzed  winds.  A  combination  of  both  may  be  desirable, 
the  former  in  middle  latitudes  and  the  latter  in  the  tropics. 

In  the  tropics  where  the  observational  errors  in  the  pressure  field  may 
be  comparable  to  the  pressure  variations  associated  with  synoptic  disturbances 
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(except  for  tropical  storms),  it  is  generally  unwise  to  attempt  to  obtain 
the  stream  function  from  the  pressure  field  by  solving  the  balance  equation, 
which,  in  fact,  is  singular  at  the  equator.   A  recommended  procedure  is  to 
calculate  the  vorticity  directly  from  the  observed  wind  field  and  then  obtain 
the  stream  function  by  solving  the  Poisson  equation  v  ij;  =  Q    .   Finally,  the 
geopotential  field  is  calculated  from  the  stream  function  with  some  form  of 
the  balance  equation.   Saha  and  Suryanarayana  (1971)  made  a  series  of  calcu- 
lations of  the  geopotential  in  this  manner  from  the  quasi-geostrophic  rela- 
tion, the  linear  balance  equation,  the  balance  equation  and  the  so-called 
vorticity  form  of  the  balance  equation  which  are,  respectively, 

V  (|)  =  fv  tj; 

V2^  =  V  •  (fViJf) 

(10) 

*v  2j<!>!)+v-  <**> 

Observed  winds  were  used  to  evaluate  the  vorticity,  Q   =  k  •  VxV  ,  and  i|f  was 
obtained  from  v  i|/  =  J  .  The  geopotential  fields  obtained  from  the  last 
three  forms  were  very  similar  and  compared  favorably  with  the  analyzed  geo- 
potential fields  at  the  850,  700,  500  and  300-rab  levels.  However,  the  last 
equation,  the  vorticity  form,  gave  the  least  RMS  error. 

b.  To  supplement  the  rotational  wind  component,  the  divergent  wind  com- 
ponent predicted  from  a  previous  analysis  may  be  used.  The  cost  in  computer 
time  amounts  to  the  solution  of  a  Poisson  equation,  V^X  =  V  •  V  ,  followed 
by  the  calculation  V  =  Vy  . 

c.  Next  apply  some  dynamic  balancing  such  as  Temperaton's  technique. 
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Step  a  and  perhaps  b  may  be  replaced  by  a  variational  method  of  analysis  in- 
cluding some  dynamical  constraints. 

Between  synoptic  times,  insert  new  data  where  available  at  intervals  of 
not  less  than  three  hours  with  an  objective  analysis  in  a  limited  region  of 
the  reports  including  first  static  balancing  to  adjust  the  wind  or  the  mass 
fields  partially  to  one  other.  This  may  be  followed  by  dynamic  balancing  to 
reduce  further  the  gravitational  noise  if  the  differencing  scheme  and  fric- 
tional  terms  are  insufficient  for  this  purpose. 

4.   Integration  Methods 

a.   Semi- imp lie  it  schemes 

Although  the  momentum  or  primitive  equations  are  simpler  and  involve 
fewer  approximations  than  the  filtered  equations,  the  presence  of  gravity 
waves  requires  a  much  smaller  time  step  to  avoid  computational  instability 
with  explicit  integration  schemes.  Otherwise  the  small  step  is  of  little 
advantage  or  even  harmful  insofar  as  meteorological  waves  are  concerned,  and 
it  is  certainly  expensive  in  terms  of  computer  time.  As  a  consequence,  there 
has  been  considerable  effort  to  circumvent  the  stability  requirement.   In  the 
Soviet  Union,  Marchuk  (1965)  introduced  a  differencing  scheme  which  treated 
the  gravity  modes  implicitly  and  the  low  frequency  meteorological  modes 
explicitly,  thus  permitting  a  much  larger  time  step.  Kwizak  and  Robert  (1971) 
successfully  applied  a  semi- implicit  differencing  method  similar  to  one 
suggested  by  Kurihara  (1965)  to  a  barotropic  500-mb  forecast.  To  illustrate 
the  technique  we  adopt  the  notation  of  Elvins  and  Sundstrom  (1973)  applied 
to  the  shallow  water  equations,  which  in  differential  form  are 
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u   =  _  (I)   -uu   -vu   +  fv 
t     Tx      x      y 

v  =  -  (b  -uv  -  v  v  -fu  (11) 

t     Ty      x      y 

(I)  =  -  (!)  (u  +  v  )  -  (uj)  +  v<b  ) 
Yt        x    y      rx    Ty 

The  pressure  gradient  terms  in  the  momentum  equations  and  the  divergence  term 
in  the  continuity  equation,  which  together  primarily  govern  the  propagation 
of  gravity  waves,  are  evaluated  semi- implicitly  while  the  remaining  terms 
are  evaluated  explicitly. 

In  the  difference  equations  u   "  and  v     in  the  continuity  equation 

are  replaced  by  substitution  from  two  momentum  equations  obtaining  an  elliptic 

An+1 
equation  in  <P 

[1-At2  f  (Dqx2  +  DQx2)]  (f+1  -  f"1)  =  F(x,y,tn,tn-1)      (12) 

Here  F  is  composed  of  terms  at  times  t  and  t   1  .   After  solving  this  equa- 
tion for   <Jr    by  one  of  the  usual  methods,  e.g.,  relaxation,  u     and 

n4-1 

v     may  be  calculated  directly  from  the  momentum  equations.   As  a  consequence 
of  the  implicit  character  of  the  difference  equations  with  respect  to  gravity 
waves,  a  much  larger  time  step  is  permitted  without  encountering  linear  com- 
putational instability.   Of  course,  extra  time  is  taken  at  each  step  to 
solve  the  elliptic  equation,  but  overall,  computer  time  is  reduced  by  a  factor 
of  3  to  4. 

Elvius  and  Sundstrom  (1973)  suggested  an  efficient  differencing 
system  which  is  staggered  in  both  space  and  time.  This  scheme  not  only 
permits  a  larger  time  step  but  also  reduces  the  phase  speed  error  of  the 
low  frequency  or  "meteorological"  mode.   They  also  developed  suitable 
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boundary  conditions  for  use  with  a  fine  mesh  model,  which  is  undergoing  fur- 
ther tests  with  realistic  initial  conditions.   Treating  the  coriolis  terms 
implicitly  permits  a  slightly  larger  time  step. 

Baroclinic  models  are  more  complicated,  but  the  semi- implicit  technique 
is  applicable  in  a  similar  manner.   Gerrity,  McPherson  and  Scolnik  (1973) 
have  developed  the  semi- implicit  difference  equations  using  Shuman's  semi- 
momentum  differencing  technique  for  the  NMC  6-layer  primitive  equation  model. 
The  model  has  run  stably  for  four  days;  however,  it  is  not  yet  operational. 

When  the  implicit  scheme  is  used,  the  phase  velocities  of  gravitational 
oscillations  are  greatly  reduced,  hence  the  geostrophic  adjustment  process 
can  be  retarded,  as  will  any  initialization  procedure  for  damping  spurious 
inertial-gravity  waves.   Experiments  by  McPherson  and  Kistler  (1973)  verified 
the  delayed  damping  of  the  gravity  waves;  however,  because  of  the  larger 
time  step,  dynamic  initialization  can  be  accomplished  with  less  computer 
time.  The  net  result  was  still  a  gain  of  2  to  1  over  the  explicit  integra- 
tion during  an  initialization  procedure. 

5.   Direct  Methods  for  Helmholtz  and  Poisson  Equations 

The  need  to  solve  Helmholtz-type  equations  in  connection  with  the  semi- 
implicit  methods  has  stimulated  interest  in  the  more  recent  "direct"  methods 
of  solving  Poisson  and  Helmholtz  equations.   Leslie  and  McAvaney  (1973)  have 
compared  the  speed  and  accuracy  of  a  number  of  methods  for  solving  equations 
of  the  form: 

vH  -  cc(x,y)  <)>  =  f(x,y)  (13) 

where  the  Helmholtz  coefficient  a  and  the  forcing  function  f  are  known. 
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When  approximated  by  the  usual  finite  difference  analogues,  the  foregoing 
equation  is  expressible  as  a  system  of  linear  equations  with  a  block  tri- 
diagonal  matrix  of  coefficients.   In  the  past,  iterative  methods  such  as 
the  Liebmann  successive  over-relaxation  method  (SOR)  have  been  widely  used 
because  of  their  simplicity.   However,  faster  direct  methods  developed  in 
recent  years  are  replacing  the  relaxation  methods,  at  least  when  a  rectangular 
region  is  involved  with  simple  boundary  conditions.   The  direct  methods  may 
be  divided  into  four  categories  as  follows: 

a.  Block  Method  which  uses  the  fact  that  A   is  block  tridiagonal. 

b.  Cyclic  Reduction  Method  (DCR)  which  reduces  the  dimensions  of  the 
matrix  to  be  solved  in  a  recursive  manner.   This  method  is  restricted 
to  certain  numbers  of  interior  points,  such  that  with  an  MXN  grid 
either  M  or  N  must  equal  2  -  1  . 

c.  Matrix  Reduction  which  through  coordinate  transformation  reduces 
the  problem  to  a  simpler  tridiagonal  form  that  has  an  easily  com- 
puted solution.   The  dimension  reduction  method  (DRM)  is  relatively 
simple  when  a  fast  Fourier  transform  is  available. 

d.  Finally,  in  certain  circumstances  the  fast  Fourier  transform  can  be 
applied  directly. 

Table  1  shows  some  comparative  solution  times  and  accuracies  (RE)  of 
the  DCR,  DRM  and  SOR  methods  applied  to  Poisson  and  Helmholtz  equations  on 
a  65  x  65  grid. 
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TABLE  1 
(from  Leslie  and  McAvaney) 


DCR 

DRM 

SOR 

RE 

TIME 

RE 

TIME 

RE 

TIME 

io-6 

0.89 

io"6 

1.25 

io"6 

42.3 

POLS  SON 

io-13 

0.90 

IO"13 

1.27 

io"10 

75.9 

5.3 

7.5 

IO"2 

io'4 

6.4 
11.0 

HELMHOLTZ 

io'6 

16.3 

It  is  clear  that  with  respect  to  Poisson  type  equations,  the  direct  methods 
are  far  superior.  To  solve  the  Helmholtz  equations  by  direct  methods,  the 
Helmholtz  coefficient  has  to  be  a  constant.   On  the  other  hand,  the  itera- 
tive (SOR)  methods  permit  a  variable  Helmholtz  coefficient  and  give  rapid 
convergence  when  that  coefficient  is  large.  This  technique  can  then  compare 

favorably  with  the  direct  methods,  particularly  if  a  little  less  accuracy  is 

-4 
acceptable,  say,  RE  <~  10   .  Also  the  SOR  methods,  which  are  most  readily 

generalized  to  irregular  domains  and  mixed  boundary  conditions,  are  simple 
to  program  and  require  minimal  storage  in  comparison  to  direct  methods. 
But  it  is  hoped  that  the  lack  of  flexibility  in  the  direct  methods  will  be 
overcome  in  view  of  the  importance  of  solving  Helmholtz  equations  in  connec- 
tion with  semi- implicit  methods  which  can  be  real  time  savers  with  respect 
to  computing. 

6.  Global  Grids 

As  computing  capability  has  improved  in  the  last  several  years  efforts 
have  increased  toward  operational  global  forecasting  and  also  in  fine  mesh 
models.  The  absence  of  lateral  boundaries  in  a  global  model  is  an  important 
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advantage  since  such  unnatural  barriers  cause  errors  in  hemispheric  models 
that  eventually  propagate  from  the  fictitious  boundaries  in  the  tropics  into 
middle  latitudes.   Numerical  forecasting  in  the  tropics,  which  admittedly 
has  severe  limitations  at  present,  would  be  quite  hopeless  if  artificial 
boundary  conditions  are  imposed  in  low  latitudes.   On  the  other  hand,  a 
tropical  band  with  middle  latitude  walls  is  not  satisfactory  either  because 
these  boundaries  are  far  too  active.   So  the  only  alternative  appears  to  be 
global  models  despite  the  vast  areas  of  little  or  no  conventional  data  in 
the  Southern  Hemisphere,  albeit  the  weather  satellites  are  helping  to  over- 
come the  data  problem. 

A  significant  difficulty  with  global  models  is  the  lack  of  a  suitable 
plane  projection  which  does  not  seriously  distort  some  areas.  The  most 
natural  approach  is  a  latitude- longitude  grid;  however,  the  convergence  of 
the  meridians  poses  a  knotty  problem,  for  as  the  distance  between  equal  longi- 
tude spacing  shrinks  toward  the  poles,  a  shorter  time  step  is  needed  to  main- 
tain linear  computational  stability.  A  time  step  short  enough  to  maintain 
stability  in  polar  regions  is  exceedingly  wasteful  in  low  latitudes.  Various 
techniques  have  been  used  to  overcome  this  difficulty  with  varying  degrees 
of  success.  The  most  common  approach  of  late  has  been  to  filter  out  or 
dampen  the  short  waves  that  would  lead  to  instability  near  the  poles  so 
that  a  relatively  large  time  step  can  be  used  throughout,  rather  than  simply 
decreasing  the  time  step  with  increasing  latitude.   In  the  Arakawa  (1972)  - 
Mintz  model,  the  procedure  consists  of  temporarily  Fourier  analyzing  fields 
which  will  be  differenced  with  respect  to  longitude  and  then  modifying  such 
Fourier  amplitude  so  that  the  CFL  stability  criterion  can  be  satisfied 
without  shortening  the  time  step.  This  principle  can  be  illustrated  with  a 
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simple  advection  equation  for  which  the  stability  criterion  for  wave  number 
k  and  maximum  wave  speed  c   is  typically  of  the  form 

^  sin  kd  <  e  (14) 

where  d.  =  a  cos  cp.  M   is  the  grid  distance  at  latitude  <P.  and  e  is  a 
constant  perhaps  unity.   By  reducing  the  amplitude  of  each  wave  component 
of  the  longitudinal  gradient  at  each  latitude  by  a  factor  S.,   the  criterion 
becomes 

<sm  >  *\   CAt  m       sin  kd  <  £  (15) 

jk  aAA  cos  cp.  v 

It  is  clear  that  for  a  fixed  &\     and  At  ,  computational  stability  can  be 
maintained  by  decreasing  S.,   as  the  latitude  and  wave  number  increase. 
This  type  of  procedure  is  applied  to  all  longitudinal  derivatives  in  the 
terms  involving  gravity  wave  propagation. 

Vanderman  (1970)  used  running  averages  of  the  tendencies  to  filter  high 
frequencies  components;  however,  this  can  cause  computational  instability. 
At  the  NOAA  Geophysical  Fluid  Dynamics  Laboratory  Holloway,  Spelman  and 
Manabe  (1973)  applied  space  filtering  to  all  time  integrated  variables  at 
each  time  step.  The  filtering  limits  the  east-west  wavelength  at  all  lati- 
tudes to  the  distance  of  two  gridlengths  at  the  equator.  The  minimum  wave- 
length is  given  by  L  .   =  4na  N   ,  where  N  is  the  number  of  gridpoints 

mxn 

around  a  latitude  circle  and  a  is  the  radius  of  the  earth.  The  maximum 
number  of  waves  at  latitude  0  is 


2TT*  COS  9  =  (N/2)  cos  6  (16) 

min 
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which  is  rounded  to  the  nearest  integer.  This  is  accomplished  by  a  Fourier 
analysis  followed  by  synthesis  with  only  the  desired  waves  at  each  latitude 
circle  included.  The  procedure  has  no  significant  effect  on  quadratic  con- 
servation properties  of  the  differencing  schemes.  The  components  of  vector 
variables  are  first  transformed  into  polar  stereographic  coordinates  to 
avoid  the  problems  associated  with  averaging  vectors  from  widely  separated 
longitudes  where  the  unit  vectors  differ  substantially  in  direction  as  dis- 
cussed by  Shuman  (1970). 

When  tested  on  barotropic  and  baroclinic  models  the  foregoing  procedure 
proved  to  be  superior  to  the  Kurihara  grid  which  has  a  poleward  decrease  in 
the  number  of  gridpoints  per  latitude  circle  in  such  a  way  as  to  keep  the 
distance  between  gridpoints  from  decreasing  appreciably.  The  problem  of 
spurious  anticyclogenesis  at  the  poles  associated  with  a  Kurihara  grid  was 
diminished.   It  was  concluded  that  fast  Fourier  methods  and  new  parallel  com- 
puters can  provide  the  necessary  speed  to  handle  the  extra  gridpoints  in 
polar  regions  when  the  simpler  differencing  schemes  are  used  with  spherical 
coordinates . 

Some  tests  were  also  conducted  by  Williamson  and  Browning  (1973)  with 
global  grids  which  verified  this  conclusion.  They  found  with  a  grid  that 
is  uniform  in  a  curvilinear  coordinate  system,  the  accuracy  of  approxima- 
tions involving  curvilinear  velocities  is  less  near  the  singularities.   In 
order  to  avoid  the  small  time  step  associated  with  a  uniform  grid,  they 
tried  the  method  of  skipping  points  near  the  poles  to  maintain  a  more  nearly 
uniform  distance  between  gridpoints,  but  the  skipped  grid  resulted  in  large 
errors.  More  accurate  interpolations  did  not  help  this  matter.  However, 
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by  applying  the  Fourier  technique  to  remove  short  wavelengths  the  errors 
were  comparable  to  a  uniform  grid  requiring  a  much  shorter  time  step. 

7.   Fourth-Order  Differencing 

In  another  aspect  of  computation- time  vs  accuracy,  it  has  been  shown 
that  a  greater  reduction  in  phase  error  can  be  achieved  per  unit  of  extra 
computer  time  by  selective  use  of  fourth-order  space  differencing  than  by 
reducing  the  grid  size  [see,  for  example,  Williams  (1972)].  This  may  be 
illustrated  with  the  simple  equation  (4).   Suppose  At  is  chosen  to  main- 
tain linear  computational  stability  in  a  P.E.  model  which  permits  a  fast 
external  gravity  wave  with  a  phase  speed  C   of  perhaps  300  m/sec  or  more 
and  C  At/Ax  °*  1  .   For  the  slower  meteorological  waves,  with  phase  speed, 
say  C   ,  the  ratio  C  At /Ax  will  be  much  less  than  one,  perhaps  a  tenth, 
or  so.   Some  computations  were  made  for  a  centered  difference  form  of  (4) 
with  both  second-  and  fourth-order  space  differences  and  second-order  time 
differencing  (leapfrog).  As  an  illustration  of  the  improvement  in  phase 
speed  accuracy,  assume  C  At/Ax  =  0.2  .  Then  the  ratios  of  the  finite 
difference  phase  speeds  to  the  true  speed  for  several  wavelengths  are  as 
follows: 

L  4Ax  6Ax  8Ax  10Ax  12Ax 
2nd  order  0.64  0.84  0.91  0.94  0.96 
4th  order    0.86   0.97   0.99   1.00    1.00 

Although  this  is  a  much  simplified  illustration,  improvements  in  phase 
speed  accuracy  of  5  to  20$  or  more  can  be  expected  with  fourth-order  space 
differences.  Moreover,  the  latter  need  only  be  applied  to  the  terms  govern- 
ing the  propagation  of  meteorological  waves  and  not  the  terms  involving 
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gravity  wave  propagation.   It  should  be  mentioned  that  the  fourth-order 
approximation  required  a  somewhat  smaller  time  step,  perhaps  20  to  25$, 
to  maintain  linear  computational  stability  [see  Haltiner  and  Williams  (1973) 
On  the  other  hand,  halving  the  grid  distance  in  the  simple  one-dimensional 
model  described  above  would  quadruple  the  computation  time,  yet  the  re- 
sulting improvement  in  phase  speeds  would  be  roughly  the  same  as  going 
from  second  to  fourth  order,  as  may  be  seen  by  comparing  the  4 Ax  and  8 Ax 
ratios  above  (0.64  and  0.91),  or  the  6 Ax  and  12 Ax  values,  0.84  and  0.96. 

8.  Mesh  Models 

Because  of  the  enormous  range  of  scales  of  atmospheric  phenomena  and 
limitations  of  even  the  most  modern  computers,  it  is  obviously  impossible 
to  model  numerically  all  phenomena  on  a  single  mesh  of  uniformly  spaced 
gridpoints.   Scales  which  are  too  small  to  be  represented  with  a  specified 
length  must  be  parameterized  in  terms  of  the  large-scale  variables  if 
their  influence  is  to  be  felt.  On  the  other  hand,  influences  from  outside 
the  region  of  integration  can  be  imposed  through  the  boundary  conditions, 
but  both  procedures  may  be  inadequate  at  times.   It  is  noteworthy  that 
some  important  mesoscale  phenomena  may  be  infrequent  in  occurrence  and  be 
predictable  for  only  a  short  range  of  time.  Consequently,  it  is  desirable 
to  superimpose  a  fine  mesh  grid  (FMG)  on  a  coarse  mesh  (CMG)  covering  a 
much  larger  area,  perhaps  a  hemisphere  or  the  entire  globe.  Quite  a  few 
meteorological  organizations,  both  foreign  and  domestic,  are  currently 
carrying  out  numerical  integrations  on  fine-mesh,  limited  area  grids. 

One  of  the  most  critical  problems  in  dealing  with  limited  area  fore- 
casts, including  the  superposition  of  different  grids,  is  the  treatment 
of  the  boundary  conditions.  This  problem  is  not  really  new  and  had  to  be 


30 


faced  in  the  first  numerical  prediction  experiments  by  Charney,  Fjortoft 
and  von  Neumann  (1950).  They  concluded  correctly  that  the  optimal  pro- 
cedure was  to  specify  precisely  as  many  boundary  values  as  required  by 
the  corresponding  linearized  equation  to  have  a  well-behaved  solution. 
Additional  values  needed  for  the  finite  difference  equations  should  be 
computed  by  extrapolations  from  interior  values.  Apparently  their  method 
of  extrapolation  proved  to  be  unstable,  as  later  shown  by  Platzman  (1954), 
and  they  simply  specified  all  values  on  the  boundary.  This  gave  stable 
results  which,  however,  were  less  accurate  and  more  stringent  than  neces- 
sary.  Having  to  maintain  constant  values  along  inflow  boundaries  can  lead 
to  large  errors  propagating  into  the  forecast  region;  however,  this  is 
certainly  less  of  a  problem  when  fine-mesh  boundaries  are  permitted  to 
change  periodically  through  the  use  of  coarse-mesh  forecasts.  Neverthe- 
less, the  latter  situation  is  in  a  sense  a  special  case  of  the  limited 
area  forecast,  for  although  the  fine-mesh  boundary  values  are  no  longer 
constant,  it  is  necessary  to  obtain  them  by  interpolation  in  space  and 
time  from  coarse-grid  values.  The  objective  then  is  to  do  this  in  such 
a  way  as  to  avoid  any  instability  at  the  boundaries  and  to  obtain  the  most 
accurate  forecast  possible. 

The  proper  procedure  for  a  barotropic  primitive  equation  model, 
proposed  by  Elvius  and  Sundstrom  (1973),  following  Charney  (1962),  is  to 
prescribe  at  all  boundary  points  the  quantity  V  -  cp£/cp  ,  which  corres- 
ponds to  the  outgoing  characteristic.  This  permits  gravity  waves  to  leave 
the  region  rather  than  be  reflected.   In  addition,  the  tangential  compo- 
nent of  velocity  is  prescribed  at  inflow  points.  Over-specification  of 
boundary  values  is  less  accurate  and  may  lead  to  parasitic  waves  and 
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perhaps  instabilities.   For  a  baroclinic  system  the  situation  is  much  more 
complicated.   Sundstrom  (1973)  recommends  an  approximation  which  specifies 
the  tangential  velocity  at  inflow  points  as  above  while  at  outflow  points 
the  tangential  velocity  is  extrapolated  from  the  interior.   Next  a  com- 
bination of  normal  velocity  component  and  potential  temperature  corres- 
ponding to  the  inward  external  gravity  wave  is  specified  on  all  boundaries, 
and  the  combination  corresponding  to  the  first  outward  gravity  wave  is 
computed  by  extrapolation  from  interior  values  in  such  a  way  as  to  avoid 
instability.  This  would  not  allow  the  internal  gravity  waves  to  leave  the 
region. 

Returning  now  to  actual  experiments  with  nesting  of  grids,  the  early 
efforts  of  superimposing  a  fine  mesh  on  a  coarse  mesh  consisted  of  inte- 
grating the  coarse  mesh  model,  for  say  24  hours,  and  saving  the  "history 
tape"  of  data  at  every  time  step.  Then  the  boundary  values  for  the  fine 
mesh  model  are  interpolated  in  space  and  time  from  the  CMG  predictions, 
which  is  an  overspecification  and  less  accurate  but  nevertheless  is  compu- 
tationally stable  [Elvius  and  Sundstrom  (1973)].  Here  the  FMG  predictions 
are  influenced  by  the  large  scales  depicted  on  the  CMG,  but  not  vice  versa. 
Nevertheless,  the  fine  mesh  gives  better  resolution  of  physical  features 
and  permits  smaller  scale  disturbances  to  develop.   Some  examples  of  this 
approach  are  the  experiments  of  Hill  (1968),  Wang  and  Halpern  (1970). 

A  successful  meteorological  experiment  which  permitted  interaction 
between  the  CMG  and  FMG  through  simultaneous  integrations  was  carried  out 
by  Harrison  and  Elsberry  (1972).   By  utilizing  a  movable  fine  mesh  centered 
over  a  traveling  one-dimensional  gravity  wave  the  FMG  produced  forecasts 
comparable  to  those  obtained  by  a  fine  mesh  everywhere.  The  boundaries 
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of  the  adjacent  CMG  and  FMG  overlapped  as  illustrated  below. 


I   I   I 


2     3    4  5   6  7  8   9  10     11     12    13 


In  their  scheme  the  FMG  integration  determines  the  values  at  points  6  and  8; 
then  these  values  are  used  for  subsequent  predictions  for  points  4  and  10  in 
the  CMG.   Similarly,  the  FMG  predictions  for  points  5  and  9  utilize  the 
values  at  points  4  and  10  determined  by  the  CMG  integrations.  Thus,  for 
example,  with  a  simple  advective  equation  bu/bt  +   U  ^u/^x  =  0  ,  the  inter- 
face equations  are 

c*u. 

ST  =  "  2d  (u6  "  u3>  •  and  (17) 

sr  -  -  t  (ue  -  V  (18) 

The  computed  changes  for  the  larger  time  steps  at  the  CMG  points  which  form 
the  boundaries  of  the  FMG  (i.e.  points  4  and  10)  are  proportioned  equally 
to  supply  the  intermediate  temporal  values  of  u,  which  are  needed  at  the 
smaller  time  intervals  to  integrate  the  adjacent  FMG  points  (i.e.  points  5 
and  9).   Similarly,  linear  spatial  interpolation  is  performed  to  obtain 
values  between  CMG  points  as  needed  for  the  FMG  integration.  As  a  conse- 
quence of  this  procedure  the  FMG  integration  influenced  the  CMG,  and  very 
importantly  the  integration  procedure  is  stable  and  created  no  significant 
noise  or  near-discontinuities  at  boundaries. 
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Applying  this  procedure  to  an  idealized  two-dimensional  tropical  disturb- 
ance and  keeping  the  fine  mesh  grid  centered  over  the  disturbance  gave  results 
equivalent  to  a  fine  mesh  everywhere  in  terms  of  forecast  intensity  and  also 
for  energy  fluxes  across  the  boundaries  dividing  the  grids.   However,  some 
"noise"  did  develop  at  the  boundary.   This  was  suppressed  by  smoothing  across 
the  boundary  and  the  addition  of  a  small  diffusion  term. 

Phillips  and  Shukla  (1973)  considered  the  two  strategies  of  one-way  inter- 
action using  the  history  tape  and  the  two-way  interaction  such  as  the  one  just 
described.  By  a  heuristic  argument  based  on  the  method  of  characteristics, 
they  inferred  that  the  two-way  interaction  procedure  would  give  a  more  faith- 
ful reproduction  of  the  proper  transmission  of  information  into  and  out  of 
the  fine-mesh  region.   Some  numerical  tests  with  the  shallow-water  equations 
showed  that  the  two-way  strategy  did  indeed  lead  to  less  error.  They  also 
found  that  the  Lax  Wendroff  two-step  scheme  gave  a  somewhat  larger  error  at 
12  hours  than  did  the  leapfrog  scheme,  but  the  reverse  was  true  at  24  and  48 
hours . 

Ookochi  (1972)  combined  a  fine  mesh  with  a  coarse  mesh  for  the  integration 
of  barotropic  primitive  equations  in  flux  form  on  a  staggered  grid.  The  re- 
sults were  essentially  a  composite  of  complete  fine-mesh  and  coarse-mesh 
integrations  with  no  significant  noise  at  the  boundaries.  The  principal 
integral  properties  involving  mass,  total  energy,  etc.  were  well  conserved 
during  the  96 -hour  experiment. 

Harrison  (1973)  describes  some  further  experiments  with  systems  of  two 
and  three  nested  grids  for  the  simulation  of  a  tropical  storm  by  integration 
of  the  primitive  equation  on  a  four-level  model.  His  calculations  demonstrated 
the  feasibility  and  advantages  of  nested  grids  in  savings  of  computer  time. 
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Presumably  better  phase  speeds  would  be  achieved  as  a  consequence  of  the 
fine  mesh  because  of  smaller  truncation  errors.  Note,  however,  that  any 
wave,  particularly  those  that  are  poorly  represented  on  the  coarse  grid, 
will  change  phase  speed  when  passing  through  the  interface  into  the  fine 
mesh,  and  again  later  when  it  leaves  the  fine  mesh.  As  a  consequence, 
erroneous  interaction  can  occur  with  that  part  that  remains  in  the  coarse 
net. 

9.  Vertical  Coordinates 

The  most  commonly  used  vertical  coordinate  in  primitive  equation  pre- 
diction models  is  the  sigma  coordinate  (a)  which  leads  to  a   -  1  and 
da/dt  =  0  on  the  lower  boundary.  This  simplification  of  the  lower  boundary 
conditions  is  accompanied  by  a  more  complicated  expression  for  the  pressure 
gradient  force. 

The  use  of  potential  temperature  as  a  vertical  coordinate  [Eliassen 
(1962)]  has  received  new  attention  in  the  last  few  years.   If  there  is  no 
heating  the  isentropic  surfaces  will  move  as  material  surfaces.  This  feature 
is  very  helpful  for  resolving  frontal  zones  and  sharp  upper  level  jets. 
A  possible  disadvantage  of  the  isentropic  coordinates  is  the  fact  that  the 
potential  temperature  shows  a  large  variation  along  the  lower  boundary. 
Primitive  equation  integrations  have  been  carried  out  by  Eliassen  and 
Raustein  (1967),  Gall  (1972)  and  Shapiro  (1973)  with  isentropic  coordinates. 
They  encountered  no  major  difficulties,  and  they  obtained  good  simulations 
of  the  jet  stream  structure.   Gall  (1972)  included  latent  heat  in  his  study. 
Bleck  (1973  a,b)  made  forecasts  with  real  data  with  the  quasi-geos trophic 
form  of  the  isentropic  potential  vorticity  equation.  The  forecasts  showed 


35 


skill  in  predicting  cyclone  development.   The  isentropic  models  show  con- 
siderable promise  for  limited  area  models  which  treat  the  smaller  scale 
synoptic  features.   These  models  have  not  been  tested  in  global  formulations 
or  for  longer  period  forecasts. 

10.  Forecasting  Skill 

With  the  greater  sophistication  in  numerical  modelling  and  the  large 
advances  in  computer  technology,  it  may  be  rightly  asked  whether  there  have 
been  concomitant  improvements  in  forecasting  skill.  To  answer  this  question 
we  may  turn  to  verification  statistics  published  by  the  National  Weather 
Service  which  are  probably  indicative  of  other  groups  as  well.   Figure  1 
from  Shuman  (1972)  shows  the  S..  scores,  which  are  approximately  a  measure 
of  the  normalized  RMS  vector  error  of  pressure  gradient,  for  30-hour,  sea- 
level  (upper  scale)  and  500-mb  (lower  scale)  forecasts  for  North  America 
from  1948  to  1971.   (Shuman  states  that  in  terms  of  practical  skill,  scores 
of  0.30  at  sea  level  and  0.20  at  500  mb  are  nearly  perfect  while  scores  of 
0.80  at  sea  level  and  0.70  at  500  mb  are  considered  essentially  worthless.) 
There  is  clearly  a  general  downward  trend  over  the  years  indicating  increas- 
ing skill.  The  improvement  in  the  latter  half  of  the  50 's  and  early  60 's 
is  ascribed  to  the  introduction  and  continuing  improvements  in  the  barotropic 
numerical  500-mb  forecasts,  while  at  sea  level  there  was  increasing  skill 
on  the  part  of  prognosticators  in  the  use  of  the  500-mb  forecasts. 

The  first  successful  baroclinic  model,  a  three-level  filtered  version, 
became  operational  in  1962,  which  is  not  reflected  clearly  in  the  graphs 
because  the  years  64-66  proved  to  be  a  more  difficult  period  to  forecast  and 
barotropic  forecasts  suffered  as  well.  Nevertheless,  a  new  plateau  of  skill 
was  established  and  there  was  a  general  downward  trend  of  S,  scores.  The 
first  surface  baroclinic  model  was  a  simplified  graphical  type  due  to 
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Dr.  Richard  Reed  which  went  into  operation  in  1962  and  drew  from  the  in- 
dependently-made numerical  500-mb  predictions.  The  continuing  improvement 
in  surface  predictions  from  1962  to  1966  was  largely  a  result  of  Reed's 
model . 

The  NMC  6-layer  primitive  equation  has  been  under  continuous  develop- 
ment since  1966;  it  is  the  first  numerical  model  at  NMC  to  produce  a  better 
sea-level  prediction  than  manual  predictions  without  NWP  guidance.   However, 
with  the  numerical  product  in  hand,  man  can  improve  the  prognosis  by  about 
five  points  on  the  average  in  terms  of  S,  scores.  The  net  result  has  been 
a  continued  improvement  in  the  skill  scores. 

Figure  2  from  Cooley  and  Derovin  (1972)  shows  the  NMC  average  36-hour, 
500-mb  S,  scores.  The  light  shading  shows  the  human  improvement.   Figure  3 
shows  the  30-hour  surface  scores. 

Considering  the  marked  increase  in  accuracy  of  surface  and  500-mb  prog- 
noses the  corresponding  improvement  in  forecasts  of  weather  elements  is 
somewhat  disappointing.   Precipitation-no  precipitation  forecasts  for  the 
periods  12-24  and  24-36  hours  improved  only  about  5$  between  1959  and  1970. 
There  was  somewhat  greater  improvement  in  temperature  forecasts.   For  ex- 
ample, the  annual  number  of  maximum  temperature  forecast  errors  equal  to 
10°F  or  greater  at  Salt  Lake  City  decreased  from  60  in  1949  to  22  in  1971. 

Sanders  (1973)  recently  reported  on  six  years  of  forecasting  tempera- 
ture and  precipitation  by  staff  and  students  of  Massachusetts  Institute 
of  Technology  for  the  local  observation  site  at  Logan  Airport.   Despite 
continuous  improvement  in  predicted  synoptic  patterns  at  the  surface  and 
aloft,  there  was  no  increase  of  skill  in  temperature  and  precipitation 
forecasts  as  measured  by  the  incremental  accuracy  over  the  climatological 


37 


mean.   As  may  be  expected,  skill  decreased  with  increasing  length  of  fore- 
cast and  more  rapidly  with  precipitation  than  with  temperature.  Minimum 
skill  occurred  in  summer,  probably  due  to  the  greater  influence  of  meso- 
scale  phenomena.   Sanders  also  indicated  that  the  lack  of  improvement  in 
forecasting  temperature  and  precipitation  and  even  a  slight  downward 
trend  in  the  latter  is  also  shared  by  NMC  forecasters.   He  suggests  that 
the  reason  for  the  first  day  forecasting  difficulties  may  lie  in  the  fact 
that  the  benefits  of  short-range  synoptic-scale  forecasts  of  the  mass 
field  have  been  maximized  and  that  the  source  of  the  errors  in  precipita- 
tion and  temperature  forecasts  are  primarily  a  result  of  mesoscale 
phenomena  such  as  fronts,  land-sea  circulations,  convection,  urban  in- 
fluences, etc.  According  to  Sanders'  results  skill  in  precipitation  fore- 
casts drops  to  10$  in  2.5  days  and  in  four  days  for  temperature.  He 
suggests  that  something  of  a  breakthrough  in  synoptic  forecasting  is 
needed  to  improve  significantly  prediction  beyond  the  first  day. 

Despite  some  variations  in  the  trends  of  forecasting  skills,  it  is  safe 
to  conclude  that,  overall,  forecasting  ability  has  shown  definite  improve- 
ment since  the  advent  of  numerical  weather  forecasts.  Moreover,  further 
improvements  can  be  expected  in  the  latter  half  of  the  70 's  as  the  current 
research  efforts  in  numerical  techniques,  simulation  of  the  physical  pro- 
cesses, and  initialization  techniques  become  operational,  along  with  better 
satellite  data  and  smaller  mesh  sizes. 

We  may  now  ask  the  question  is  there  no  limit  to  the  ultimate  range 
and  accuracy  of  weather  forecasts  if  one  is  willing  to  spend  enough  money 
to  provide  the  necessary  data  and  enough  computing  power  for  the  calcula- 
tions? To  answer  this  question  the  predictability  of  the  atmosphere  must 
be  considered. 
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11.   Atmospheric  Predictability 

Lorenz  (1965)  pointed  out  that  perfect  prediction  can  never  be  expected 
because  (a)  the  governing  laws  are  not  perfectly  known  but  only  approxima- 
tions, (b)  nor  is  the  system  of  equations  strictly  deterministic  and  finally 
(c)  the  initial  state  is  not  perfectly  known.  Moreover,  he  shows  that  even 
for  a  determinate  system  of  equations,  such  as  one  normally  used  in  numeri- 
cal weather  prediction,  separate  solutions  which  are  nearly  identical  at 
the  initial  time  do  not  remain  nearly  identical  as  time  progresses,  but 
eventually  become  as  different  as  two  solutions  chosen  randomly  at  the  same 
time  and  day  of  the  year.   The  atmosphere  contains  some  periodic  oscilla- 
tions -  principally  the  annual  and  diurnal  variations  and  their  overtones 
-  which  are  predictable  at  essentially  infinite  range.   However,  accurate 
long  range  prediction  of  the  remaining  features  is  impossible  because  the 
initial  state  of  the  atmosphere  will  always  be  imperfectly  known. 

Lorenz  (1969a)  suggested  three  approaches  that  might  be  used  to  esti- 
mate the  predictability  of  the  atmosphere: 

a.  The  dynamical  approach  wherein  a  system  of  equations  closely  resem- 
bling the  atmospheric  equations  is  integrated  from  slightly  different 
initial  conditions  and  then  the  rate  of  amplification  of  the  differences 

is  determined; 

b.  The  empirical  approach  which  utilizes  similar  weather  types,  usually 
referred  to  as  "analogues"  and  determines  their  subsequent  separation  in 
statistical  terms  as  a  function  of  time; 

c.  The  empirical-dynamical  approach  which  utilizes  a  new  set  of  equa- 
tions (derived  from  the  atmospheric  equations)  which  describes  a  spectral 
distribution  of  forecast  errors. 
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Most  of  the  studies  which  follow  (a)  used  general  circulation  models 
which  had  been  integrated  until  they  were  in  approximate  balance.   At  this 
time  another  integration  was  begun  with  an  initial  state  which  consisted 
of  the  balanced  field  plus  a  small  departure.  As  this  solution  departed 
from  the  control  experiment  the  growth  rate  of  the  error  field  and  the 
range  of  predictability  were  obtained.  Experiments  by  Mintz  and  Arakawa, 
Smagorinsky  and  Leith  have  been  summarized  by  Charney,  et  al  (1966). 
Charney  estimated  the  RMS  doubling  time  of  small  errors  to  be  about  4  to 
5  days  and  a  predictability  limit  imposed  by  typical  observational  errors 
of  about  two  weeks.   Smagorinsky  (1969)  presented  the  results  of  experi- 
ments carried  out  at  the  NOAA  Geophysical  Fluid  Dynamics  Laboratory  with 
their  10- level  primitive  equation  general  circulation  model.  A  random 
temperature  disturbance  with  a  standard  deviation  of  0.5  C  at  all  levels 
was  added  to  the  control  run.  The  vertical  average  of  the  standard  devia- 
tion dropped  from  0.5  C  to  0.2  C  after  one  day,  reflecting  a  "geostrophic 
adjustment"  between  the  initial  disturbed  wind  field  and  the  undisturbed 
wind  field.  Thereafter,  the  error  growth  was  exponential  for  the  next 
seven  days  with  a  doubling  time  of  about  2.5  days.   Smagorinsky  concluded 
that  the  deterministic  limit  of  predictability  for  synoptic  scale  disturb- 
ances is  about  three  weeks;  however,  the  current  practical  limit  is  about 
one  week.  He  also  noted  that  short-wave  predictability  decays  most  rapidly. 

The  cause  of  the  exponential  growth  of  small  errors  in  these  studies 
was  attributed  to  baroclinic  instability  [Charney,  et  al  (1966)].  However, 
Lorenz  (1972)  suggested  that  the  error  growth  might  be  due  to  barotropic 
instability  of  wave  disturbances.   In  particular  he  investigated  the 
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stability  of  a  finite  amplitude  unbounded  Rossby  wave.   He  found  instability 
when  the  waves  are  sufficiently  strong  or  the  wave  number  is  sufficiently 
high,  and  the  growth  rates  are  comparable  to  the  growth  rates  obtained  from 
predictability  studies.   Lilly  (1973)  emphasized  the  importance  of  baro- 
tropic  instability  in  a  predictability  study  which  employed  a  2-dimensional 
model.  Recently,  Lorenz  (1973b)  has  concluded  that  baroclinic  instability 
is  the  most  important  cause  of  lack  of  predictability  in  the  atmosphere. 

Lorenz  (1969c)  utilized  approach  b  to  obtain  the  rate  of  separation  of 
two  fields  which  were  initially  similar.   Five  years  of  twice-daily  height 
values  of  the  200- ,  500- ,  and  850-mb  surfaces  on  a  grid  of  1003  points  were 
obtained.  A  weighted  root-mean-square  height  difference  is  used  as  a  measure 
of  the  difference  between  two  states,  or  the  error.   For  each  pair  of  states 
occurring  within  one  month  of  the  same  time  of  year,  but  in  different  years, 
the  error  was  computed.  Numerous  mediocre  analogues  were  found  but  there 
were  no  really  good  ones.  The  smallest  errors  had  a  doubling  time  of  about 
eight  days.   Since  larger  errors  grow  less  rapidly  this  is  probably  an  over- 
estimate of  the  doubling  time.   Extrapolation  with  the  aid  of  a  quadratic 
hypothesis  indicates  that  very  small  errors  double  in  about  2.5  days.  This 
compares  very  well  with  the  numerical  results  reported  by  Smagorinsky  (1969). 

In  a  more  recent  paper  Lorenz  (1973a)  has  used  the  same  data  set  to  in- 
vestigate the  range  of  predictability.   States  of  the  atmosphere  separated 
by  12  days  or  less  are  found  on  the  average  to  resemble  each  other  more 
closely  than  randomly  selected  states,  even  after  an  adjustment  for  the 
seasonal  trend  has  been  made.   Higher  correlations  were  obtained  with  a 
form  of  damped  persistence.  These  results  demonstrate  the  existence  of 
partial  predictability  of  instantaneous  weather  patterns  at  least  12  days 
in  advance. 
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Lorenz  (1969b)  followed  approach  c  with  a  statistical  treatment  of  the 
2 -dimensional  vorticity  equation.   He  derived  an  equation  for  the  "error 
energy"  which  is  obtained  from  the  difference  between  two  solutions. 
Thompson  (1957)  considered  similar  equations.   The  linearized  equations 
are  written  in  spectral  form  and  ensemble  averages  are  taken.  A  further 
assumption  is  included  to  close  the  set  of  predictive  equations.   The  equa- 
tions require  the  specification  of  a  mean  energy  spectrum,  and  Lorenz  used 
the  -5/3  law  for  the  smaller  scales.  The  equations  were  integrated  numeri- 
cally from  an  initial  error  distribution;  the  error  energy  of  each  wave 
number  was  not  allowed  to  exceed  the  mean  atmospheric  energy  for  that 
wave  number.  The  numerical  solutions  show  a  rapid  growth  of  the  error  for 
the  very  small  scales  and  a  very  slow  growth  for  large  scales.  Lorenz 's 
solutions  show  that  the  range  of  predictability  for  cumulus  scales  is 
almost  an  hour,  while  the  synoptic-scale  motions  can  be  predicted  a  few 
days  ahead.   Predictability  for  the  largest  scale  disappears  after  about 
17  days.   In  particular  if  the  initial  error  is  confined  to  the  smallest 
scales  the  error  in  those  scales  rapidly  reaches  the  maximum  error.  Then 
growth  occurs  in  the  next  larger  scales  until  the  error  is  propagated  to 
the  largest  scale.   Lorenz  points  out  that  if  the  error  in  the  smallest 
scale  is  reduced  then  the  range  of  predictability  of  the  largest  scale 
features  will  only  be  increased  by  a  time  interval  which  is  less  than  the 
range  of  predictability  of  the  small  scales  where  the  error  was  reduced. 
Although  these  results  are  dependent  on  the  closure  assumption  and  the 
energy  spectrum  which  is  used  the  study  shows  clearly  the  importance  of 
the  nonlinear  propagation  of  errors  between  different  scales  of  motion. 

Robinson  (1967,  1971)  assumed  that  the  dynamic  equations  do  not  allow 
one  to  predict  motions  of  a  given  scale  over  periods  longer  than  the  fluid 
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elements  of  this  scale  maintain  their  identity  against  turbulent  diffusion 
by  smaller  scale  motions.   Then  based  upon  the  dissipation  characteristics 
of  the  atmosphere,  he  derived  predictability  ranges  of  a  few  days  for 
synoptic-scale  motions  and  about  one  hour  for  cumulus  scale  motions,  which 
are  roughly  the  same  as  obtained  by  Lorenz  (1969b). 

Leith  (1971)  and  Leith  and  Kraichman  (1972)  have  also  used  homogeneous 
isotropic  turbulence  models  to  study  atmospheric  predictability.  Leith 
(1971)  considered  two  dimensions  and  closed  his  equations  with  the  eddy- 
damped  Markovian  approximation.   Leith  and  Kraichman  (1972)  considered  both 
two  and  three  dimensions  and  they  used  Kraichman's  (1971)  closure  method 
which  is  based  on  the  test-field  model.   The  solutions  in  three  dimensions 
showed  that  even  where  there  is  a  strong  energy  cascade  there  is  also  an 
upscale  propagation  of  error.  The  two-dimensional  solutions  used  a  "-3" 
power  energy  spectrum  for  the  smaller  scales.   Charney  (1971,  1973)  showed 
that  quasi-geos trophic  three-dimensional  flow  should  have  a  "-3"  power 
energy  spectrum  under  proper  conditions.  Leith  and  Kraichman  (1972)  found 
that  the  two-dimensional  flow  was  more  predictable  than  three-dimensional 
flow.  They  concluded  that  an  initial  state  determined  with  a  horizontal 
resolution  feasible  with  a  satellite-based  observing  system  would  result 
in  significant  predictability  of  large  scale  motion  for  more  than  one  week. 
They  also  point  out  that  the  test-field  model  probably  underestimates  rather 
than  overestimates  predictability  times. 

In  summary,  the  studies  of  predictability  suggest  that  specific  weather 
patterns  and  events  are  predictable,  or  partially  so,  only  for  a  period  of 
at  most  several  weeks;  however,  the  possibility  of  predicting  general  trends 
of  perhaps  temperature  and  precipitation  above  or  below  normal  for  longer 
periods  is  not  precluded  as  yet. 
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12.  Stochastic  Dynamic  Prediction 

Recent  studies  in  stochastic  dynamic  prediction  are  closely  related  to 
atmospheric  predictability.   Epstein  (1969)  developed  a  method  for  includ- 
ing the  effect  of  random  errors  in  initial  conditions  in  a  forecast  model. 
The  forecast  equations  which  he  treated  were  of  the  spectral  type  with 
quadratic  nonlinear it ies .  He  first  took  an  ensemble  average  of  the  basic 
equations  which  leads  to  predictive  equations  for  the  ensemble  average  or 
the  expected  value  of  each  dependent  variable.   These  equations  also  in- 
clude the  variances  and  covariances  between  the  dependent  variables.   Equa- 
tions for  the  time  change  of  the  latter  quantities  which  are  second  moments 
are  obtained  by  multiplying  each  equation  by  each  of  the  variables  and 
carrying  out  the  ensemble  averages.  These  equations,  however,  contain 
the  third  moments.   Equations  for  the  time  change  of  the  third  moments  are 
obtained  in  the  same  way  and  they  contain  the  fourth  moments.  This  is 
merely  the  closure  problem  of  classical  turbulence  theory. 

Epstein  closed  his  equations  by  neglecting  the  third  moments  about  the 
instantaneous  mean.  He  integrated  the  resulting  equations  for  the  Lorenz 
(1960)  three  component  system.   Initially  he  specified  the  expected  value 
of  each  variable  as  well  as  its  variance  due  to  possible  observational 
errors.  The  initial  covariances  were  set  equal  to  zero.  A  deterministic 
forecast  was  also  performed  with  original  prediction  equations.  Numerical 
integrations  showed  that  the  deterministic  solution  eventually  diverges 
from  the  stochastic  prediction  of  the  expected  value.   The  latter  should  be 
a  better  forecast  in  the  statistical  sense.  The  predicted  variances  grow 
in  time  which  gives  the  influence  of  the  initial  error.  Also  a  large  number 
of  deterministic  forecasts  were  made  with  slightly  different  initial  condi- 
tions. These  numerical  solutions  were  then  averaged  using  a  Monte  Carlo 
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method.   The  difference  between  the  Monte  Carlo  average  and  the  approximate 
stochastic  solution  is  a  measure  of  the  error  which  is  caused  by  the  closure 
assumption.   In  some  of  the  solutions  these  quantities  remained  close  while 
in  others  they  began  to  depart  after  a  few  days. 

Fleming  (1971a)  has  continued  this  work  with  a  two-level  quasi-geostrophic 
spectral  model  similar  to  the  one  used  by  Lorenz  (1965).   For  interpretation 
he  divided  the  energy  into  a  "certain"  energy  and  "uncertain"  energy.   He 
also  reconsidered  Epstein's  closure  assumption  which  neglected  the  third 
order  moments.   Fleming  tested  the  quasi-normal  approximation  which  computes 
the  fourth  order  moments  in  the  third  order  equations  in  terms  of  the  second 
order  moments.   This  requires  the  numerical  solution  of  the  third  order  equa- 
tions in  addition  to  the  others.  He  found  that  this  approximation  was  better 
up  to  a  certain  time,  after  which  it  became  unstable.  He  then  considered  a 
modified  form  in  which  a  linear  damping  term  was  added  to  the  third  order 
equations.  This  eddy-damped  quasi-normal  approximation  was  stable  and  gave 
excellent  results  in  most  cases. 

Fleming  (1971b)  used  the  eddy-damped  quasi-normal  approximation  in  a 
stochastic  formulation  of  the  barotropic  model.  He  used  two  wave  numbers 
in  latitude  and  15  in  longitude  in  a  study  of  predictability.   The  predicta- 
bility times  obtained  are  similar  to  those  obtained  by  Lorenz  (1969b).   Pre- 
dictability studies  using  the  baroclinic  model  of  the  earlier  paper  showed 
that  predictability  is  increased  when  heating  and  friction  are  present. 
Fleming  (1972)  considered  the  effect  of  random  variations  in  the  thermal 
forcing  and  the  effect  of  errors  resulting  from  the  improper  treatment  of 
the  smaller  scales. 
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Leith  (1973)  has  concluded  that  the  stochastic  dynamic  method  cannot 
be  used  for  the  larger  numerical  models  because  of  the  very  excessive  com- 
puter time  requirements.   He  suggests  instead  the  use  of  the  Monte  Carlo 
procedure  which  involves  a  collection  of  deterministic  forecasts.   The 
procedure  was  evaluated  with  a  simple  two-dimensional  turbulence  model  of 
large-scale  atmospheric  motions.   A  study  of  the  dependence  on  sample  size 
of  mean  square  forecast  error,  both  with  and  without  a  final  linear  regres- 
sion correct  ion, showed  a  considerable  improvement  with  moderate  sample 
sizes  over  conventional  single  forecasts  without  regression.  These  results 
suggest  that  a  sample  size  of  about  10  may  be  adequate  for  producing  the 
error  variance  information  needed  for  optimal  data  assimilation. 
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13.  Spectral  Methods 

The  spectral  method  expresses  the  spatial  variation  of  the  prediction 
fields  in  terms  of  a  series  of  orthogonal  functions.   The  coefficients  in 
the  series  are  now  the  forecast  quantities  rather  than  gridpoint  values  of 
the  original  dependent  variables.   We  will  illustrate  the  technique  with 
the  barotropic  vorticity  equation: 


^7  V^  +  |k  x  Vjf  .  vcv2^)  =  0  (19) 


where  ij;   is  the  stream  function  and  the  beta  term  has  been  neglected  for 
simplicity.   We  expand  ty      into  the  following  series: 


K^.^.t)  =  £  Ya(t)  Y  (x^)  .  (20) 
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The  functions  Y   are  orthogonal  and  normalized  so  that 
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V  Ya  dS  "  ^  •  (21) 


where  Y_    is  the  complex  conjugate  of  Yp  and 

p  P 
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i: 


a(3  a  =  p 


The  integration  in  (21)  is  carried  out  over  the  entire  forecast  region. 

The  functions  Y   satisfy  the  equation 
a 

*   Ya  +  Va  =  °  '  (22) 

where  the  eigenvalues  K   are  positive  and  increase  for  the  smaller  scale 
eig en functions.   Substitute  (20)  into  the  forecast  equation  (19),  multiply 
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by     Y        ,    integrate  over  the   entire  region  and  use  the  relation  (21)  which 


gives 


-K7  J+L    £     ly.o.p  *0  »p  -  Q   •  <23) 


the   interaction  coefficients,   which  can  be  computed  once  and   for  all,   are 
given  by 


L     „   q  =  -  Ko    \  Y  *  Ik   •   VY^  x  vY.   dS    . 


(24) 


Equation  (23)  represents  an  infinite  number  of  ordinary  differential  equa- 
tions when  all  appropriate  values  of  y     are  considered.   In  practice  the 
series  given  by  (20)  is  truncated  in  such  a  way  that  the  desired  meteorologi- 
cal features  are  properly  described.  The  equations  for  the  remaining 
coefficients  can  then  be  integrated  in  time  numerically. 

Lorenz  (1960)  treated  equation  (19)  with  this  procedure  in  cartesian 
coordinates.   In  that  case  the  eig en functions  are  products  of  sines  and 
cosines  and  the  interaction  coefficients  take  on  a  particularly  simple  form. 
He  noted  that  when  the  series  is  truncated  total  energy  will  still  be  con-  . 
served  if  the  excluded  coefficients  are  set  to  zero  in  the  interaction  sums. 
He  also  demonstrated  the  usefulness  of  very  low  order  systems  with  his  study 
of  a  3-component  set. 

Silberman  (1954)  first  applied  this  technique  in  spherical  coordinates 
although  he  kept  the  zonal  mean  flow  fixed  in  time.   Platzman  (1960)  and 
Baer  and  Platzman  (1961)  performed  a  complete  treatment  of  the  barotropic 
vorticity  equation  in  spherical  coordinates  and  also  carried  out  some 
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experiments.   In  spherical  coordinates  the  eigenf unctions  are  spherical 
harmonics  which  are  composed  of  sines  and  cosines  of  the  longitude  and 
Legendre  polynomials  of  the  sine  of  the  latitude.   With  these  functions 
the  interaction  coefficients  are  more  complicated  and  they  have  more  non- 
zero elements  than  in  cartesian  coordinates. 

The  spectral  method  has  several  advantages  over  the  gridpoint  method. 
First  the  spectral  method  computes  spatial  derivatives  exactly  so  that  the 
phase  error  which  occurs  with  the  finite  difference  method  is  eliminated. 
Also  the  aliasing  that  occurs  with  finite  differences  is  excluded  and  as 
a  result  it  is  easy  to  conserve  certain  quantities  which  are  conserved  in 
the  continuous  equations.   Poisson  equations  are  easily  solved  without 
relaxation  or  other  inversion  techniques  because  of  relation  (22).  Another 
advantage  is  the  treatment  of  global  motions  without  the  presence  of 
singularities. 

The  most  important  disadvantage  of  the  spectral  method  is  that  it 
requires  much  more  computer  time  than  the  gridpoint  method,  if  there  are 
very  many  degrees  of  freedom.  This  can  be  seen  from  equation  (23)  which 
shows  that  for  each  degree  of  freedom  many  products  must  be  computed  in 
the  nonlinear  term.   With  the  gridpoint  method,  for  each  degree  of  freedom, 
there  are  only  a  few  products  involving  quantities  at  surrounding  grid- 
points.   The  spectral  method  is  generally  more  complicated  to  apply  and 
it  is  not  convenient  for  complicated  boundaries.  The  spectral  method  also 
suffers  when  values  must  be  evaluated  locally  such  as  in  determining  con- 
densation criteria.  We  shall  see  that  some  of  these  difficulties  have 
been  alleviated  with  recent  developments. 
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Robert  (1966,  1968)  modified  the  spectral  technique  in  two  studies  with 
the  primitive  equations  on  the  sphere.   In  place  of  the  spherical  harmonic 
functions  he  used  the  following  function  of  sines  and  cosines: 


G   (A,cp)  =  e  n  cos1  '  cp  sin  c 


cp  sin  cp  .  (25) 

n 


These  functions  are  not  orthogonal,  but  the  spherical  harmonic  functions 
can  be  written  as  a  sum  as  of  the  functions  given  by  (25).  The  nonlinear 
terms  are  easily  computed  with  these  functions,  and  the  number  of  elements 
involved  is  much  less  than  with  the  spherical  harmonics.  When  the 
orthogonality  relations  are  required  the  spherical  harmonics  can  be  formed 
in  terms  of  these  functions.  Using  this  simplified  procedure  Robert 
carried  a  number  of  spherical  integrations  with  a  relatively  small  number 
of  terms. 

A  very  important  advance  for  the  spectral  method  was  the  development  of 
the  fast  Fourier  transform  by  Cooley  and  Tukey  (1965).  Their  technique 

allows  Fourier  transformation  of  a  field  with  N  degrees  of  freedom  with 

2 
order  N  log  N  operations,  while  the  direct  method  requires  order  N 

operations.  This  allows  rapid  calculation  of  local  fields  by  summing 

the  series  with  the  fast  Fourier  transformation,  such  as  might  be  required 

for  condensation  tests. 

Orszag  (1969,  1970)  has  used  the  fast  Fourier  technique  to  save  compu- 
tation time  in  the  evaluation  of  the  nonlinear  terms.  He  carries  out  all 
differentiation  of  quantities  while  they  are  represented  by  the  series, 
but  products  of  fields  are  performed  with  the  fields  at  gridpoints.  This 
requires  an  inverse  Fourier  transform  to  obtain  data  at  gridpoints  and  a 
regular  Fourier  transformation  to  return  to  spectral  form.  The  fast  Fourier 
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transformation  can  be  used  in  both  of  these  operations  with  a  great  time 
savings.   When  this  procedure  can  be  used  it  greatly  reduces  prediction 
time,  for  large  numbers  of  variables,  because  the  sum  in  (23)  is  replaced 
by  a  much  faster  operation. 

Merilees  (1973a)  has  developed  an  algorithm  for  computing  the  sum  of 
a  series  of  spherical  harmonics.  The  algorithm  is  approximately  10-20 
times  faster  than  the  standard  method,  but  it  suffers  a  precision  problem 
and  it  breaks  down  when  the  resolution  is  too  high.   When  this  method  is 
used. following  Orszag,  for  rapid  computation  of  the  nonlinear  terms,  the 
prediction  time  for  a  global  spectral  model  can  be  greatly  reduced. 

Bourke  (1972)  formulated  a  global  spectral  model  based  upon  the  one- 
layer  shallow  water  equations.   He  used  the  vorticity  and  divergence  equa- 
tions instead  of  the  2  components  of  the  equation  of  motion.   The  technique 
for  evaluating  the  nonlinear  terms  which  was  developed  by  Orszag  (1970) 
and  Eliasen,  et  al  (1970)  is  employed.   The  computational  efficiency  of 
the  model  was  found  to  be  far  superior  to  that  of  an  equivalent  model  based 
on  the  interaction  coefficients.  This  model,  in  integrations  of  116  days, 
satisfied  the  principles  of  conservation  of  energy,  angular  momentum  and 
square  potential  vorticity  to  a  high  degree. 

Daley  (1973)  used  Bourke1 s  (1972)  model  to  examine  the  possibility  of 
using  different  spatial  resolution  for  different  dependent  variables.  He 
suggested  that  for  the  smaller  scales  it  would  be  sufficient  to  predict 
only  the  wind  field  because  for  those  scales  the  pressure  field  adjusts  to 
the  wind  field.   In  a  test  the  smaller  scales  were  predicted  with  a  filtered 
model  and  the  larger  scales  with  the  full  primitive  equations.   He  found 
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that  this  combined  system  gave  better  forecasts  than  a  primitive  equation 
system  with  intermediate  resolution.   The  application  of  this  procedure  to 
baroclinic  equations  may  be  more  difficult. 

Recent  developments  have  made  the  spectral  models  much  more  competi- 
tive with  gridpoint  models  with  respect  to  computational  efficiency.   For 
the  same  resolution  the  spectral  models  can  be  expected  to  give  better 
forecasts,  and  the  forecasts  may  be  better  even  with  somewhat  poorer  reso- 
lution. A  series  of  tests  are  needed  to  determine  whether  or  not  better 
forecasts  can  be  made  with  the  spectral  method  with  the  same  computational 
effort. 

Orszag  (1972)  and  Merilees  (1973b)  have  proposed  a  forecast  technique 
which  is  a  combination  of  the  finite  difference  method  and  spectral  method. 
This  technique  which  is  known  as  the  pseudospectral  approximation  employs 
finite  differences  except  in  the  computation  of  spatial  derivatives. 
Before  the  derivatives  are  computed,  the  gridpoint  data  is  fast  Fourier 
analyzed.  Continuous  derivatives  are  then  computed  and  the  series  is  summed 
with  the  fast  Fourier  transform.  This  procedure  is  used  for  both  longi- 
tudinal and  latitudinal  derivatives.  The  latitudinal  derivatives  in  some 
cases  require  a  sign  change  at  the  poles.  Merilees  (1973b)  found  excellent 
results  in  a  test  prediction  of  Haurwitz  waves  with  the  shallow  water 
equations. 

14.  Finite  Element  Method 

The  extensive  success  of  the  finite  element  method  in  solid  mechanics 
has  attracted  the  interest  of  investigators  in  other  fields.   Salinas  (1973) 
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has  given  a  brief  review  of  the  method  which  proceeds  as  follows:   At  the 
start,  an  approximate  solution  in  the  form  of  a  linear  combination  of  speci- 
fied (basis)  functions  is  assumed.   The  coefficients  (multipliers)  of  the 
basis  functions  are  to  be  determined  to  yield  the  best  approximate  solu- 
tions.  This  is  accomplished  by  minimizing  a  measure  of  the  error  (called 
the  residual  function)  associated  with  the  assumed  solution.   Several  tech- 
niques are  available  for  minimizing  the  residuals.  Normally  the  basis 
functions  are  defined  in  such  a  way  that  they  have  a  simple  variation 
(quadratic,  cubic,  etc.)  over  an  element  of  area  (piecewise  polynomials), 
outside  of  which  they  are  zero.  The  "elements"  can  have  a  variety  of 
shapes.  The  finite  element  method  is  a  generalization  of  the  method  of 
weighted  residuals. 

Newton  (1973)  has  successfully  applied  the  technique  to  gravity  waves 
radiating  out  from  an  oscillating  ship.   Gallagher  and  Chan  (1973)  have 
treated  the  steady  state  circulation  in  Lake  Ontario.  Triangular  and 
otherwise  shaped  elements  were  used  to  accurately  fit  the  lake's  coastline. 
Wang,  et  al  (1972)  applied  the  method  to  the  one-dimensional  shallow  water 
equations.  They  obtained  better  results  for  both  the  linear  and  nonlinear 
cases  than  with  the  usual  finite  difference  methods.   Price  and  MacPherson 
(1973)  have  applied  the  technique  to  the  two-level  general  circulation  model 
developed  by  Mintz  and  Arakawa  (Langlois  and  Kwok,  1969).  They  arranged 
the  elements  in  such  a  way  that  the  area  of  each  element  is  nearly  con- 
stant over  the  entire  globe.  They  also  provided  for  a  subregion  in  which 
the  elements  telescoped  to  a  smaller  size.   Predictions  with  this  model 
gave  encouraging  results. 

For  meteorology,  the  principal  advantage  of  the  method  as  shown  by 
Price  and  MacPherson  (1973),  is  the  possibility  of  easily  changing  the 
element  size  and  shape.   This  would  be  useful  in  situations  where  meshed 
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models  are  now  used  and  also  on  the  sphere  where  the  elements  could  be  kept 
at  a  fixed  scale. 

15 .  Summary 

Research  on  initialization  has  intensified  in  the  last  few  years.  Most 
operational  forecast  models  employ  the  primitive  equations  which  are  very 
sensitive  to  the  initial  balance  between  the  wind  and  the  pressure  fields. 
Short  range  prediction  of  precipitation  is  also  very  sensitive  to  the  initial 
divergence  field.  The  development  of  global  models  requires  proper  initiali- 
zation in  the  tropics  which  may  differ  from  middle  latitudes.  The  most 
promising  initialization  techniques  include  the  variational  formulations 
and  the  averaging  of  backward  and  forward  forecasts  about  the  initial  state. 
The  assimilation  of  nonsynoptic  data  into  a  forecast  model  has  received  con- 
siderable attention.  This  has  been  motivated  by  availability  of  nearly  con- 
tinuous satellite  data.   The  main  problem  here  involves  controlling  the 
imbalance  introduced  into  the  forecast  by  the  localized  new  data. 

Several  groups  have  developed  limited  area  models  which  are  designed  to 
give  better  forecasts  in  specific  regions.   These  models  use  boundary  condi- 
tions from  hemispheric  or  global  models  and  in  some  cases  the  limited  area 
forecast  is  allowed  to  affect  the  exterior  region.   Isentropic  models  have 
been  developed  and  applied  to  smaller  scale  synoptic  features. 

Operational  predictions  of  pressure  and  wind  have  shown  continuing  im- 
provement as  a  consequence  of  progressively  better  numerical  models.  How- 
ever, the  concomitant  improvement  in  forecasting  precipitation  and  tempera- 
ture has  not  kept  the  same  pace.   It  is  expected  that  forecasts  of  weather 
elements  will  improve  as  the  limited  area  models  produce  better  descriptions 
of  mesoscale  systems. 


54 


Studies  of  predictability  suggest  that  specific  weather  patterns  and 
events  are  predictable,  or  partially  so,  only  for  at  most  several  weeks. 
The  possibility  of  predicting  general  trends  above  or  below  normal  for  longer 
periods  is  not  as  yet  precluded.   Perhaps  the  consideration  of  atmosphere- 
ocean  interaction  will  lead  to  longer  range  prediction  of  weather  trends. 

The  stochastic  dynamical  prediction  procedure  is  closely  related  to 
studies  of  atmospheric  predictability.  The  procedure  uses  the  initial  un- 
certainty in  the  initial  state  to  predict  the  uncertainty  in  the  forecast 
at  later  times  as  well  as  the  expected  value.  This  method  cannot  be  used 
for  operational  forecasts  in  the  near  future  because  of  the  large  computer 
time  requirements.   However,  an  indication  of  the  growth  of  uncertainty 
can  be  obtained  by  examining  a  relatively  small  number  of  deterministic 
forecasts  with  slightly  different  initial  conditions. 

The  development  of  the  semi-implicit  method  is  a  major  advance  in  finite 
differencing.  The  method  treats  the  gravity  wave  terms  implicitly  and  the 
advection  terms  explicitly.  This  permits  the  use  of  a  longer  time  step 
with  the  additional  requirement  that  an  elliptic  equation  be  solved  at  each 
time  step.  The  net  effect  of  this  is  a  reduction  in  computer  time  by  a 
factor  of  at  least  2.  This  has  lead  to  the  development  of  new  methods  for 
solving  elliptic  equations  which  are  faster  and  more  accurate  than  the 
traditional  relaxation  methods. 

Present  operational  forecasts  often  underpredict  the  movement  of  synoptic 
disturbances  and  this  error  has  been  attributed  to  space  truncation  error 
in  the  finite  difference  equations.   Several  groups  are  experimenting  with 
fourth  order  space  differencing  in  order  to  reduce  this  error.  Also  the 
pseudospectral  method  has  been  put  forward  as  a  method  to  eliminate  space 
truncation  effects. 
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New  interest  has  developed  in  the  use  of  the  spectral  method  since  the 
introduction  of  the  fast  Fourier  transform.  The  use  of  this  transform 
greatly  speeds  up  computation  time  with  the  spectral  method.   Further 
tests  are  required  to  determine  whether  the  spectral  method  will  give 
better  forecasts  than  the  finite  difference  method  for  the  same  amount 
of  computer  time. 

The  finite  element  method  of  solid  mechanics  has  been  recently  intro- 
duced into  meteorology.  This  method  has  the  advantage  of  a  very  flexible 
element  size,  but  it  has  not  been  widely  tested  on  meteorological  problems, 
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B.  APRIL  1958  TO  FEBRUARY  1964    BASED  ON  NWP  500-MB  BAR0TR0PIC 

AND  3-LEVEL  BAR0CLINIC  GUIDANCE 
C    MARCH  1964  TO  MAY  1966    BASED  ON  REED  SURFACE  GUIDANCE. 

D.  JUNE  1966  TO  AUGUST  1967    BASED  ON  P.E.  SURFACE  GUIDANCE 

E.  SEPTEMBER  1967  TO  DECEMBER  1971.  BASED  ON  P.E.  SURFACE  GUIDANCE 

D  INDICATES  HUMAN  IMPROVEMENT. 
Figure   2.--NMC   average   30-hour   surface   S,    scores 
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